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1.0  INTRODUCTION 


1.1  BACKGROUND 


In  analyzing  and  characterizing  experimental  or  numerical  electro¬ 
magnetic  response  data,  one  generally  desires  to  extract  parameters  that 
can  be  related  to  physical  characteristics  of  the  system  being  studied. 
One  set  of  physically  related  parameters  is  the  complex  natural  reson¬ 
ances  of  the  system  and  their  related  coefficients.  Indeed  these  reson¬ 
ances  have  been  recognized  as  important  for  some  time.  Recently  the.  con¬ 
cept  of  electromagnetic  resonances  has  been  applied  to  problems  of  par¬ 
ticular  interest  to  the  Navy. 


In  1971  Baum  [1]  developed  the  formalism  known  as  the  Singularity 
Expansion  Method  (SEM)  which  enables  one  to  write  any  electromagnetic 
response  of  a  system  as  an  expansion  of  the  complex  resonances  or  poles 
and  residues  of  that  system.  Baum  developed  this  formalism  for  the 
electromagnetic  pulse  (EMP)  community  so  that  the  external  current  dis¬ 
tribution  on  aircraft  and  ground  support  systems  could  be  characterized 
concisely.  It  was  not  until  the  singularity  expansion  method  came  into 
being  that  it  was  possible  to  determine  the  modal  resonances  and  the 
excitation  coefficients  of  each  mode  for  a  structure  with  an  arbitrary 
incident  excitation. 


In  1974  Mains  and  Moffatt  [2]  introduced  the  concept  of  using  the 
complex  natural  resonances  of  a  body  as  a  basis  for  target  recognition. 
They  made  use  of  the  fact  that  a  few  natural  resonances  of  a  body  are 
adequate  to  distinguish  the  body  within  a  finite  collection  of  bodies. 
They  also  made  use  of  the  knowledge  that  the  natural  resonances  of  a 
body,  as  manifested  in  a  scattered  waveform,  are  aspect  independent.  The 
identification  procedure  which  was  used  was  to  first  obtain  a  set  of 
multi-frequency  radar  scattering  data  from  an  unknown  target.  A  pre- 
diction-correlator  type  of  processing  was  then  applied  to  select  the 
real  target  from  a  catalog  of  candidate  targets  and  their  resonances. 
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In  both  Baum's  and  Mains  and  Moffatt's  early  work  the  resonances 
were  obtained  from  a  set  of  equations  which  characterized  the  electro¬ 
magnetic  response  of  the  body,  much  as  a  circuit  theorist  finds  his 
resonances  by  solving  a  differential  equation.  Many  shapes  have  been 
numerically  analyzed  to  date  to  obtain  their  characteristic  resonances. 
However  for  complex  shapes  and  configurations,  it  is  of  great  interest  to 
be  able  to  obtain  the  resonances  from  experimental  data.  In  particular 
obtaining  the  resonance  parameters  from  transient  response  data  from  EMP 
simulators  and  transient  radar  ranges  is  of  keen  interest  to  the  military 
community. 

Spectral  electromagnetic  response  data  usually  lends  itself  to  the. 
visual  identification  of  these  natural  frequencies.  The  damping  constants 
cannot  be  obtained  as  easily,  however,  and  often  are  calculated  from  es¬ 
timates  of  the  quality  factor  (Q) .  Similarly,  temporal  response  data 
generally  allows  one  to  visually  determine  the  dominant  natural  frequency 
in  a  response.  If  enough  data  are  present,  the  damping  constant  of  this 
dominant  frequency  can  also  be  determined.  However,  temporal  data  are 
usually  Fourier  transformed  to  the  frequency  domain  so  that  the  higher 
order  modes  can  be  identified  visually.  Visual  identification  of  the 
natural  resonances  of  a  system  is  not  ideal  by  any  means,  particularly  if 
the  system  has  many  low  Q  modes.  Hence  a  numerical  resonance  extraction 
procedure  is  very  much  in  demand. 

About  seven  years  ago  Prony's  algorithm  [3]  was  applied  for  the  first 
time  to  the  problem  of  numerically  extracting  the  natural  resonances  from 
transient  electromagnetics  response  data.  The  first  application  of  Prony's 
method  was  to  the  numerically  generated  transient  current  on  a  thin  dipole. 
The  results,  which  were  reported  at  the  1974  USNC/URSI  meeting  [4]  by 
Mittra  and  Van  Blaricum  gave  a  set  of  resonances  (poles)  which  compared 
very  closely  with  the  first  ten  even  modes  previously  calculated  by 
Tesche  [5].  As  a  result  of  this  initial  demens tration,  several  researcn- 
ers  began  studying  Prony's  method  to  determine  its  utility  for  analyzing 
several  kinds  of  transient  data  and  to  look  for  solutions  to  some  of  the 
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problems  inherent  in  the  process.  In  addition,  Brittingham,  Miller  and 
Willows  [6]  demonstrated  that  a  procedure  parallel  to  the  time  domain 
Prony's  method  could  be  applied  to  frequency  domain  data. 

Some  of  the  initial  questions  which  were  asked  about  Prony's  method, 
ar ^  which  are  in  part  still  begin  studied,  were: 

1.  Will  Prony's  method  work  if  multiple  poles  are  present? 

2.  How  does  one  determine  a  priori  the  order  of  the  system? 

That  is,  how  many  poles  are  contained  in  the  response  data? 

3.  What  effects  do  noisy  data  have  on  the  Prony  algorithm? 

4.  How  do  we  insure  or  know  the  accuracy  of  the  poles  returned? 

Thca-'  questions  were  all  addressed  to  some  extent  in  Van  Blaricum's 
c-^rfertat  [ ?].  It  was  found  that  Prony's  method  would  work  for  the 

C-  ii  c.  mu.ti.ple  poles  without  any  change  in  the  pole  searching  algorithm. 
Tv.  thods  werr  discovered  by  which  the  order  of  the  system  could  be 
deter- The^o  methods  are  the  Householder  orthogonalization  proce¬ 
dure  arii.'  the  Eigenvalue  method.  Examples  of  these  methods  can  be  found 
in  Vin  Blaricum's  thesis  [7]  and  appear  in  the  special  EMP  issue  of  AP-S 
L 8 j .  A  preliminary  scudy  of  the  effects  of  noise  on  the  pole  extraction 
algorithm  and  the  order  determination  algorithms  was  presented  in  [7]  and 
[8]  also.  However  the  problem  of  noise  and  Prony's  method  is  a  very 
complex  one  which  has  not  been,  even  at  this  time,  completely  answered. 
Several  alternatives  to  Prony's  method  have  been  suggested  and  investi¬ 
gated  with  the  hope  of  finding  a  cure  for  all  noise  problems.  Among 
these  alternatives  are  variations  on  Prony's  method  [9-10];  the  pencil- 
of-functions  method  [11]  which  is  presented  in  Appendix  D;  iterative 
generalized  least  squares  presented  in  Appendix  A;  column  Prony's  method. 
Appendix  F;  the  adaptive  method,  Appendix  C;  and  Evan's  and  Fischl's 
method.  Appendix  G. 

To  understand  the  different  procedures  for  pole  extraction  and  the 
difficulties  involved  it  is  necessary  to  discuss  the  main  elements  of  the 
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existing  procedures.  The  process  of  (extracting  the  natural  resonances 
and  their  associated  residues  from  a  transient  signal  has  four  main  steps 
as  shown  in  Figure  1.1. 

The  first  step  is  the  determination  of  the  order  of  the  system.  At 
this  step  one  decides  how  many  poles  the  system  response  function  has  so 
that  the  proper  model  order  can  be  obtained.  It  has  been  found,  through 
a  combination  of  parameter  studies  and  trial  and  error,  that  if  the  order 
of  the  system  is  underestimated  then  the  extracted  poles  will  deviate  sub¬ 
stantially  from  the  true  poles.  Similarly  if  the  order  of  the  system  is 
overdetermined,  the  algorithm  produces  extraneous  poles.  The  presence 
of  the  extraneous  poles  causes  the  residues  of  the  true  poles  to  be 
inaccurate  and  also  results  in  unnecessary  computation  time.  The  presence 
of  noise  in  data  makes  the  determination  of  the  system  order  a  very  complex 
problem.  Up  to  this  time  many  methods  have  been  used  to  determine  the 
order  but  they  either  break  down  when  noise  is  present  or  they  are 
dependent  on  trial  and  error  or  the  intervention  of  the  user.  For  analysis 
of  massive  amounts  of  data,  as  In  the  case  of  EMP  data,  or  for  radar  target 
identification  a  totally  automated  method  is  a  must.  In  Section  3  (Volume 
I)  of  this  report,  a  procedure  capable  of  automatically  determining  the 
proper,  order  without  any  knowledge  of  the  noise  level  is  presented. 


Once  the  order  of  the  system  has  been  determined  the  coefficients 
of  the  linear  predictor  equation  or  Prony's  difference  equation  must  be 
solved.  The  degree  of  difficulty  of  this  step  depends  on  the  noise  level 
in  the  data  and  on  the  proper  determination  of  the  order  of  the  system. 

In  Section  4,  certain  factors  which  greatly  effect  the  accuracy  of  the 
coefficients  are  discussed. 

Once  the  difference  equation  coefficients  ate  obtained,  the  roots 
th 

of  anN  order  polynomial,  N  being  the  system  order,  must  be  found. 

Many  root  finding  routines  exist  but  Muller's  method  [13]  appears  to  be 


FIGURE  1.1.  The  Elements  of  the  Extraction  of  Natural 
.'.sonances  and  Associated  Residues  from  a 
Transient  Signal 


the  optimal  method.  While  this  is  a  key  step  in  the  procedure  it  is 
totally  dependent  or.  the  accuracy  of  the  coefficients  which  were  obtained 
in  the  pre'^ous  step. 

The  final  step  i3  the  solution  of  the  residues  which  are  associated 
with  the  system  poles  or  singularities.  These  residues  are  obtained  by 
solving  a  simple  linear  matrix  equation.  In  many  problems,  such  as 
target  identification,  the  residues  are  not  required  and  hence  this  is 
certainly  not  a  critical  step. 

1.2  RECENT  WORK  IN  RESONANCE  ESTIMATION 


Kulp  [14]  has  recently  studied  the  effect  of  sampling  rate  on  the 
accuracy  of  Prony's  method.  The  sampling  rate,  we  have  found,  is  one  of 
two  factors  that  determine  the  performance  of  Prony's  method  with  noi3y 
data.  The  other  factor  is  statistical  bias  in  the  estimates.  If  boch 
of  these  problems  are  treated,  Prony's  method  performs  admirably  well 
with  noisy  data. 

Cures  for  the  bias  problem  and  the  problem  of  extraneous  resonances 
ha^e  been  proposed  by  Henderson[15] .  The  technique  involves  the  use  of 
eigenvalue  decomposition  to  construct  the  coefficients  of  Prony's 
difference  equation.  Although  we  cannot  prove,  as  yet,  that  the  use  of 
eigenvalue  decomposition  allows  bias-free  estimation  of  the  coefficients, 
we  have  seen  indications  that  this  is  the  case.  In  Section  2,  Volume  I 
of  this  report,  indications  that  lead  us  to  this  bias-free  remark  are 
described. 

It  is  these  two  recent  works  plus  order  selection  techniques  developed 
by  ETI  that  set  the  stage  for  testing  of  the  automatic  resonance  extrac¬ 
tion  procedures  in  Phase  II  of  the  current  contract. 
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1.3  OBJECTIVES  AND  CONTRIBUTIONS  OF  THE  PRESENT  WORK 


Three  major  objectives  were  set  forth  at  the  beginning  of  the 
current  contract: 

1.  To  review  and  summarize  the  state-of-the-art  of  procedures 
for  extracting  resonances  from  transient  data  which  include 
Prony’s  method,  adaptations  of  Prony's  method  and  other 
techniques  of  system  identification  and  parameter  estimation. 

2.  To  investigate  methods  for  automatically  determining  the 
proper  order  of  a  system  represented  by  a  noisy  set  of  data 
and  the  limitations  of  the  procedures  for  doing  so. 

3.  To  investigate  methods  for  automatically  extracting  the 
proper  system  resonav cos  from  a  set  of  noisy  data  once 
the  proper  order  has  been  established. 

At  the  present  time  these  objectives  have  been  fully  satisfied.  Our 
present  understanding  of  the  resouance  extraction  problem  will  enable 
us  to  construct  and  test  practical  automatic  procedures  for  resonance 
extraction  in  the  second  phase  of  the  current  contract.  In  looking  for 
automated  methods,  several  other  important  issues  such  as  biased 
estimates,  iterative  versus  one  shot  processes,  and  pattern  recognition, 
to  name  a  few,  were  investigated.  Many  possible  techniques  were 
examined  in  a  search  for  automated  methods.  These  techniques  are 
summarized  in  appendix  form  in  Volume  II  of  this  report. 

1.4  SOME  REMOTE  SENSING  APPLICATIONS  OF  RESONANCE  EXTRACTION 

The  ideal  end  product  of  this  type  of  research  is  an  automatic  real 
time  technique  which  will  take  noisy  transient  signals  and  estimate  the 
complex  resonances  of  the  structure  and  the  errors  in  these  estimates. 

In  studying  existing  techniques  with  this  one  end  product  in  mind  it 
became  clear  to  us  that  the  method  developed  is  dependent  on  the  specific 
uses  intended. 


The  three  major  remote  sensing  applications  of  the  complex  resonances 
are  target  recognition,  target  intelligence,  and  target  camouflage  and  de¬ 
coying.  Target  recognition,  as  we  define  it,  assumes  that  previous  know¬ 
ledge  of  the  targets  in  question  is  known.  Target  recognition  is  also 
assumed  to  be  a  real-time  or  a  near  real-time  process.  The  poles  ex¬ 
tracted  from  the  measured  targec's  signature  can  be  compared  to  a  target 
pole  catalog  for  identification.  While  this  requires  work  in  building 
the  initial  library  it  makes  the  actual  pole  extraction  procedures  less 
complicated  (fortunately  since  they  have  to  be  real  time)  because  the 
true  poles  will  not  have  to  be  separated  from  the  noise  poles. 

Target  intelligence  assumes  that  we  are  trying  to  discover  the 
shape  and  type  of  vehicle  we  are  seeing  having  never  seen  it  before. 

Hence  this  target  will  not  be  in  a  lihrary  or  catalog.  This  requires 
that  the  true  system  poles  be  separable  from  the  noise  poles.  In  addition 
the  relation  between  pole  patterns  or  locations  and  the  target's  physical 
characteristics  must  be  known.  The  process,  however,  does  not  have  to 
be  real-time  and  can  usually  have  human  input. 

Target  camouflage  and  decoying  are  potential  techniques  based  on 
knowing  the  relationships  between  pole  patterns  and  physical  character¬ 
istics.  Knowing  these  relationships  a  target  can  be  given  apparent  (to 
the  radar)  new  characteristics  by  modifying  its  measured  pole  pattern. 

Hence  one  could  make  a  cruise  missile  either  apparently  disappear  by 
overdamping  its  resonant  return  or  make  it  look  like  another  vehicle  by 
appropriately  moving  the  resonances.  Target  camouflage  is  much  like 
target  intellegence  in  its  requirements.  The  main  requirement  is  develop¬ 
ing  the  knowledge  of  the  relationship  among  the  poles  and  physical  config¬ 
uration. 

Before  much  more  work  is  done  in  this  general  theoretical  area  of 
resonance  extraction  from  transient  signals  it  will  be  necessary  to  very 
accurately  define  the  desired  type  of  system  application  to  which  the 
methods  will  be  applied.  In  Section  4,  Volume  I  of  this  report,  we  de¬ 
fine  tentative  and  very  rough  forms  of  the  procedures  to  be  applied  to 


the  target  recognition  and  target  intelligence  problems.  These  rough 
forms  will  gain  further  definition  in  the  second  phase  of  the  current 
contract. 

1.5  OUTLINE  OF  THIS  REPORT 

Volume  I  presents  major  results  and  conclusions  of  the  research. 

Section  2  addresses  the  question  of  what  the  best  methods  of  estimating 
the  resonances  for  this  application  are.  Among  the  estimation  procedures 
studied  are  the  conventional  technique  of  using  an  Inhomogeneous  set  of 
equations,  eigenvalue  decomposition  tecnniques,  Jain's  method,  and  iter¬ 
ative  techniques.  Section  3  studies  procedures  for  selecting  the  proper 
order  and  eliminating  extraneous  resonances.  One  major  and  original 
contribution  by  ETI  is  presented  in  Section  3,  namely,  the  maximum  likelihood 
procedure  for  order  selection.  This  procedure  is  capable  of  automatically 
determining  the  proper  order  without  any  prior  knowledge  of  the  noise 
level.  Also,  in  Section  3,  Henderson's  procedure  for  eliminating  extran¬ 
eous  resonances  is  tested  by  numerical  example.  The  procedures  presented 
in  Section  3  allow  an  efficient,  automatic  procedure  for  order  selection 
and  resonance  estimation  to  be  constructed  by  using  eigenvalue  decomposition 
of  the  data.  Section  4  discusses  various  practical  problems  that  must  be 
addressed  before  resonance  extraction  procedures  for  a  practical  radar 
target  recognition  system  can  be  designed.  This  section  serves  to  relate 
the  concepts  in  Sections  2  and  3  to  the  application  under  consideration 
in  the  current  contract.  Section  4  addresses  the  basic  problem  of  what 
is  the  best  sampling  rate  and  period  of  observation  for  the  purposes  of 
resonance  extraction.  Possible  procedures  and  systems  for  both  target 
recognition  and  target  intelligence  using  natural  resonances  are  presented. 
Also  in  Section  4,  two  key  problems  with  Prony's  method  are  identified. 
Methods  for  effectively  treating  these  problems  are  presented  which  allow 
Prony's  method  to  be  effective  at  high  noise  levels. 
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Volume  II  consists  of  appendices  which  present  the  details  of 
several  techniques  that  have  been  investigated  during  this  contract. 
Volume  I  makes  references  to  Volume  II  for  the  details  of  specific 
procedures. 

Volume  III  contains  a  large  bibliography  of  Prony*s  method.  Also 
in  Volume  III  is  a  translation  of  Prony's  original  paper. 


2.0  COMPARISON  OF  PROCEDURES  FI/R  RESONANCE  ESTIMATION 


2.1  DEFINITION  OF  THE  RESONANCE  EXTRACTION  PROBLEM 

In  this  report  we  consider  the  problem  of  estimating  the  resonances 
of  a  linear,  time- Invariant,  physical  system,  such  as  an  electromagnetic 
scatterer,  from  the  measured  time-domain  system  response  to  the  measured 
time-domain  excitation.  We  assume  that  a  distributed  physical  system 
such  as  a  scatterer  can  be  adequately  modeled  by  a  single-input,  single¬ 
output  linear  system  where  the  input  may  be  the  incident  electric  field 
time  history  and  the  output  is  the  reflected  electric  far-field  time 
history. 

In  general,  both  the  excitation  and  response  are  required  to  estimate 
the  resonances  of  the  system.  However,  for  the  first  part  of  this  section 
we  consider  the  case  where  the  excitation  is  a  delta  function  so  that  the 
impluse  response  can  be  observed.  In  practical  situations  the  impulse 
response  can  usually  be  obtained  by  deconvolution.  (Appendix  H  presents 
one  deconvolution  technique.)  The  measured  impulse  response  can  be 
expressed  as 


A 

“  wi  +  ei  "  2j  Aj  exP  <8j  +  ei  - 


J-l 


j-1 


A  z  +  e 
J  i 


(1) 


for  i«0,l...,N  where  z.  -  exp(s.At)  are  the  resonances  or  poles  in  terms 

8  J  J 


of  z-transform  variable,  the  a,  are  the  poles  in  terms  of  the  Laplace- 

J  th 

transform  variable,  the  Aj  are  termed  "residues",  qi  denotes  the  i  sample, 


Ng  is  the  number  of  samples,  and  At  is  the  time  increment  between  sucessive 


samples.  Equal  increments  are  assumed. 
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Distributed  systems  normally  have  an  infinite  number  of  modes,  i.e., 
k»oo .  In  practice,  the  number  of  poles  that  are  modeled,  n,  need  not  be 
infinite  but  only  sufficiently  large  to  yield  an  adequate  approximation  to 
the  measured  response.  The  "best"  value  of  n  is  strongly  dependent  on  the 
number  of  dominant  resonances  and  the  noise  level. 


Many  procedures  are  available  for  estimating  the  resonances.  But, 
for  the  first  part  of  this  section,  we  consider  only  the  class  of  procedures 
based  on  Prony's  difference  equation  which  has  the  form 

n 

2  a  w  *0,  i-0,l,...,m.  ^2) 

J-Q  j  j+i 

where  denotes  the  unocrrupted  value  of  If  is  replaced  by  q1  the 

equation  becomes 
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where  is,  in  general,  nonzero,  and  is  termed  the  "equation  error'  .  In 
matrix  form  (2)  becomes  Wx“0  and  (o)  becomes  Qx—d 
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are  the  (M  x  H)-dimensional  data  matrices, x  -  ^a0*ai » *  *  '  is  the 

dimensional  parameter  vector,  and  d  -  tdu’d1 . drJT  is  the  M-dimensional 

equation  error  vector.  For  convenience  we  use  M-m+1  and  N-n+1 . 
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The  central  problem  of  all  noniterative  resonance  extraction  procedures 
is  choosing  the  parameter  vector  to  minimize  the  magnitude  of  the  equation 
error.  Various  procedures  for  minimizing  d  are  available.  In  this  section, 
we  examine  some  of  these  procedures  and  isolate  one  procedure  that  is  the 
most  appropriate  for  resonance  extraction. 

Once  x  has  been  adjusted  to  minimize  d,  the  roots  z  of  the  equation 


E  ^  ^  ■  o  (4) 

j-0  2 


are  estimates  of  the  resonances  in  the  z-plane  and  define  estimates  of  the 


s-plane  resonances  through  s^ 


£nz*/At,  for  j n  where,  for  the 


moment,  n«k  is  assumed.  Estimates  of  the  residues  can  be  taken  as  the 
coefficients  that  minimize  the  "true  error"  defined  a3 


m+n  n 

E  [q.  “  £  a;  (z-)1]2  .  (5) 

i-0  1  j-0  2  2 


The  minimizing  A'  can  be  chosen  by  using  a  standard  least-squares  technique. 


An  important  observation  can  be  made  from  (4):  If  the  ctj  are  all 
multiplied  by  the  same  constant,  the  roots,  zj*  ,  are  unchanged*  From  this 
observation  it  is  concluded  that  the  magnitude  of  the  x  vector  is  irrelevant 
in  estimating  the  resonances.  Only  the  direction  of  x  is  important. 


In  this  section  we  assume  that  the  true  number  of  poles, k,  is  finite 
and  known  and  that  n-k.  The  section  on  order  selection  procedures  discusses 
what  courses  of  action  are  available  when  k  is  not  known  and  possibly  infinite 
The  presence  of  measurement  error  produces  a  data  matrix  of  rank  N  while  the 
true  rank  without  measurement  error  is  k.  Measurement  error,  then,  greatly 
complicates  estimation  of  the  true  order. 
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2.2  QUADRATIC  ERROR  MINIMIZATION 


One  procedure  for  minimizing  the  equation  error  is  straight  minimization 

of  its  norm, ||d||,  or,  equivalently,  the  minimization  of  the  square  of  the 
T  T 

norm,  d  d.  We  refer  to  d  d  as  quadratic  error  (QE) .  Since  d*Qx,  the. 

trivial  solution  x*“0  will  always  provide  the  minimum  QE.  Other  values  of 

x  can  provide  this  same  minimum  value  only  if  the  rank  of  Q  is  less  than  N 

which  occurs  only  in  the  absence  of  measurement  error.  Nevertheless,  we  may 

proceed  in  the  usual  fashion  to  derive  a  set  of  normal  equations  by  equating 

T 

the  derivative  of  d  d  with  respect  to  x  to  zero: 

(dTd)  -  2QTQx  -  0 


T 

which  implies  that  the  value  of  x  that  satisfies  Q  Q  x-0  can  furnish  a  minimum. 

T 

Unfortunately,  when  measurement  error  is  present  Q  Q  is  nons.ingular,  and 
hence,  only  the  trivial  solution,  x*0,  exists  as  we  noted  above.  The  trivial 
solution  is  not  desirable  since  the  zero  vector  contains  none  of  the 
directional  information  that  is  necessary  to  estimate  the  resonances. 

One  way  out  of  this  dilemma  is  to  fix  one  element  of  the  x  vector  at  a 

T 

nonzero  value  and  adjust  the  remaining  elements  to  minimize  d  d.  Suppose 

we  let  a  -1.  Then  Qx“d  can  be  written  Qx+q-d  where  x  =  [u,,»a1,...a  .  ]T  , 

n  01  n-i 


T 

and  q-lq  ,q  , ,...q  ]T .  The  deriviative  of  d  d  with  respect  to  x  is 
n  n+l  Tn+n 


7x  (dT^)  *  V-  [(Qx  +  q)T(Qx  +  q)] 

-  2QTQx  +  2QTq  . 


Equating  the  derlviatlve  to  zero  implies  that 


-  — -1 

XT  -  -  [Q  Q]  Q  q  . 


With  or  without  measurement  error  QXQ  is  always  nonsingular,  and  hence,  its 
inverse  exists.  It  follows  that  can  always  be  used  to  find  a  nontrivial 

parameter  vector.  jL  is  referred  to  as  the  reduced  or  inhomogeneous  solution 

1  _T__  _T_ 
because  it  satisfies  the  inhomogeneous  equation:  Q  Qx^«Q  q. 

The  use  of  the  inhomogeneous  solution  has  traditionally  been  the 
procedure  for  obtaining  a  parameter  vector.  We  will  examine  this  solution 
procedure  more  closely  after  another  procedure  is  introduced. 

2.3  MINIMIZATION  OF  NORMALIZED  QUADRATIC  ERROR 

T  T 

We  have  seen  that  d  d  becomes  zero  when  x  x  is  zero  which  leads  to  the 

problem  of  the  trivial  solution.  An  error  criterion  that  does  not  necessarily 

tend  to  zero  when  the  parameter  vector  tends  to  zero  is  the  normalized 

T  T  T  T 

quadratic  error  (NQE)  defined  as  d  d/x  x.  The  deriviative  of  d  d/x  x  with 
respect  to  x  is 


x  T 

X  X 


A 


-  2x 


<*V 


Equating  the  deriviative  to  zero  implies  that 
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This  is  the  well-known  form  of  an  eigenvalue  problem.  The  only  values  of 

T 

x  capable  of  furnishing  a  minimum  of  NQE  are  the  eigenvectors  of  Q  Q.  NQE 
assumes  the  minimum  as  the  lowest  eigenvalue  when  x  assumes  the  value  of 
the  weakest  eigenvector.  We  conclude  that  minimizing  NQE  effectively  avoids 
the  trivial  solution  and  is  equivalent  to  using  an  eigenvalue  technique  to 
find  a  minimizing  parameter  vector.  The  weakest  eigenvector  may  be  used  to 
estimate  the  resonances.  Price  [9]  arrived  at  this  result  using  different 
but  equivalent  arguments. 

2.4  SINGULAR  VALUE /EIGENVALUE  DECOMPOSITION 

T 

Another  way  of  representing  the  eigenvalue  analysis  of  Q  Q  is  ithe 

T 

eigenvalue  decomposition  (EVD)  of  Q  Q  which  can  be  expressed  as 

T  T 

•  QlQ  -  VAjjV1 

where  V  is  an  (N  x  N) -dimensional  orthogonal  matrix  whose  columns  are  the 
T 

eigenvectors  of  Q  Q  and  A  is  an  (N  x  N) -dimensional  diagonal  matrix  whose 
diagonal  elements  are  the  eigenvalues  of  Q  Q. 

T 

In  the  same  way,  the  eigenvalue  analysis  of  QQ  can  be  expressed  in 
the  EVD  form: 

T  T 

QQ1  -  UAgU1 

where  U  is  an  (M  x  M) -dimensional  orthogonal  matrix  whose  columns  are  the 

T 

eigenvectors  of  QQ  and  A^  is  an  (M  x  M) -dimensional  diagonal  matrix  whose 
elements  are  the  eigenvalues  of  QQ^. 

A  more  general  form  of  EVD  for  nonsquare  matrices,  known  as  singular 
value  decomposition  [16]  (SVD) ,  effectively  performs  both  EVD’s  mentioned 
above.  The  SVD  of  Q  is  expressed  as 
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Q  -  USV" 

where  U  and  V  are  the  matrixes  that  were  previously  defined  and  S  is  an 
(M  x  N)-dimensional  matrix  whose  diagonal  elements  are  non-negative  and 

are  called  the  singular  values  of  Q.  The  singular  values  are  related  to 

T  T 

the  eigenvalues  by  ■  S  S  and  ■  SS  . 

2.5  USE  OF  EVD/SVD  TO  APPROXIMATE  THE  DATA  MATRIX 

T  T 

We  have  seen  that  straight  minimization  of  d  d  leads  to  Q  Q»0  which 

for  the  noisy  case  has  only  the  trivial  solution  since  Q  Q  is  nonsingular. 

Two  methods  to  avoid  the  trivial  solution,  fixing  an  element  of  the  parameter 

vector  and  EVD/SVD  analysis,  have  been  mentioned.  Another  interpretation 

of  EVD/SVD  analysis  promotes  understanding  of  what  these  methods  accomplish. 

EVD/SVD  analysis  can  be  interpreted  as  a  method  to  approximate  Q  Q  or  Q 

with  a  matrix  of  lower  rank  or,  in  other  words,  a  singular  matrix.  When 
T 

Q  Q  is  replaced  by  a  singular  approximant  matrix,  the  homogeneous  equation 

mentioned  above  has  a  nontrivial  solution.  It  is  shown  that  this  nontrivial 

T 

solution  is  equivalent  to  the  weakest  eigenvector  of  Q  Q. 

T 

A  singular  approximant  of  P"Q  Q  can  be  constructed  by  using  the  EVD 
of  QTQ  as 


where  A  _  is  the  diagonal  matrix  constructed  from  by  forcing  the  smallest 
diagonal  element  or  eigenvalue  to  zero.  The  eigenvector  corresponding  to 
the  smallest  nonzero  eigenvalue  must  satisfy  P  x«0  and  hence,  it  is  the 
nontrivial  solution  to  the  best  approximant  of  the  homogeneous  equation 
q^qx»0.  P  is  the  optimal  approximant  to  P  in  the  sense  that  is  the  singular 
matrix,  Ps,  that  minimizes  the  euclidean  matrix  norm  of  Pfl-P. 
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An  identical  procedure  is  the  use  of  the  SVD  of  Q  to  construct  an 
approximant  of  Q: 

ft  =*  us  vT 

s  s 

where  Sg  is  constructed  from  S  by  forcing  the  smallest  singular  value  to 
zero.  Since  £  "Q  the  resulting  solution  is  identical  to  the  solution 

S  S  S  rp 

obtained  with  EVD  of  Q  Q. 

One  important  observation  should  be  made  about  the  approximation  Qg : 
This  matrix  is  not  a  Hankel  matrix  in  general.  Therefore,  we  cannot  use 
this  matrix  to  form  an  approximation  to  the  original  waveform. 


2.6  SUPERIORITY  OF  EVD/SVD  TECHNIQUES  FOR  RESONANCE  ESTIMATION 

So  far  two  methods  to  find  a  minimizing  parameter  vector  have  been 
found: 

1.  The  inhomogeneous  solution 

2.  EVD/SVD  analysis 


In  the  following  we  show  that  the  directions  (not  just  the  magnitudes)  of 
the  parameter  vectors  obtained  using  these  two  procedures  are  different, 
and  hence,  the  estimates  of  the  resonances  differ.  We  also  show  that  in 
one  sense  EVD/SVD  analysis  is  superior  to  the  inhomogeneous  solution. 


Let  XjjQg  be  the  weakest  eigenvector  of  QiQ  which  has  been  scaled  so 
that  <*n“l-  (Note  that  scaling  the  eigenvector  does  not  change  its  direction, 
and  hence,  the  estimates  of  the  resonances  are  unchanged.)  Let  x^g  denote 
the  reduced  parameter  vector  which  is  formed  from  Then  QXjjqg  * 

^*NQE  +  q  "  dNQE‘  By  taking  the  pseudoinverse  of  Q  we  solve  for  x^^,: 
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*NQE  "  -[Q  Q]  Q  q  +  [Q  Q1  Q  V  . 

_  _ip_  _2-x- 

Since  the  inhomogeneous  solution  is  defined  as  ■  -[Q  Q]  Q  q  ,  the 
difference  between  the  two  estimates  is 


*NQE  xi  *  IQ  QJ  Q  dNQE  # 

When  noise  is  present,  this  difference  can  be  shown  to  be,  in  general, 
nonzero.  Since  an*=l  for  both  procedures,  the  parameter  vectors  constructed 
from  these  reduced  vectors  must  differ  in  direction  in  the  general  case, 

Now  that  the  two  procedures  have  been  shown  to  produce  different 
resonance  estimates,  we  provide  an  argument  to  show  that  Xjjqg  is  the 
better  of  the  two  estimates  because  it  is  perturbed  less  than  the  inhomo¬ 
geneous  solution  by  the  presence  of  noise  in  the  data. 

If  we  assume  that  the  noise  corrputing  the  waveform  is  zero-mean, 

2 

Gaussian-distributed,  and  uncorrelated  with  variance  a  then  we  can  show 
that  [7] 

$(QTQ)  “  WTW  +  IMa2  , 

5 (QTQ)  =  WTW  +  IMa2  ’ 

_T—  -T- 

and  ?(Q  q)  ■=  W  w  where  I  denote  the  identity  matrix  of  appropriate  dimension 
in  each  case,  W  and  w  are  the  uncorrupted  matrices  corresponding  to  Q 
and  q,  respectively,  £  denotes  the  expectation  operator,  and  M  is  the 
number  of  rows  of  Q  and  Q. 
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In  order  to  perform  a  completely  deterministic  analysis  of  the  effects 
of  noise  on  each  procedure,  we  replace  QTQ,  QTQ  and  QTq  with  their  expect¬ 
ations  we  derived  above. 


The  EVD  of  WTW  is  denoted  by  WTW  -  .  Then  the  EVD  of 


5[Q  Q]  is 


£[QTQ]  -  VW[A*  +  IMa2]VWT 


T  T 

Since  W  W  and  5(Q  Q)  possess  the  same  eigenvectors,  we  conclude  that 

T 

noise  has  no  net  effect,  on  the  average,  on  the  eigenvectors  of  Q  Q  including 
Xjjqj,  .  Proof  is  not  yet  available  to  support  the  claim  that  the  expectation 
of  the  eigenvectors  is  uneffected  by  noise.  However,  we  provide  the  following 
geometrical  argument:  The  eigenvectors  represent  the  extremes  of  the  directions 
present  in  the  data,  that  is,  either  the  strongest  or  weakest  directions  are 
represented  by  eigenvectors.  A  purely  random  noise  component  has  no  dominant 
directions  and  tends  to  perturb  the  eigenvectors  equally  in  all  directions 
neighboring  their  unperturbed  or  noise-free  directions.  Because  of  the  symmetry 
in  the  perturbations,  the  average  or  expectation  of  the  perturbed  eigenvectors 
should  be  equal  to  the  unperturbed  eigenvectors  themselves. 


Let  xIE>  which  is  referred  to  as  the  expected  Inhomogeneous  solution, 
satisfy 


£ (QTQ)  xte  -  5(QTq) 


or  (VTW  +  I  Mu2)  L  •  WTw  . 
n  XE 


(6) 


We  assume  that  the  expected  solution  is  representative  of  actual  solutions, 
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The  following  observations  can  then  be  made: 


1.  When  0«o>  i«e.,  no  noise,  the  expected  solution  reduces  to  the 
noise-free  or  exact  solution. 


3.  From  1  and  2,  the  expected  solution  is  highly  dependent  on 
the  noise  level. 

The  third  observation  contrasts  the  effects  of  noise  on  the  Inhomogeneous 
solution  and  on  xNqE* 

From  the  foregoing  analysis  we  conclude  the  following: 


x.,-_  and  x_  lead  to  different  resonance  estimates. 
NQE  I 


2.  Resonance  estimates  from  x^  are  considered  to  be  inferior  to 
those  from  x^pE  because  is  expected  to  display  a  strong 
dependence  on  the  noise  level  while  is  expected  to  remain 

invariant  with  the  noise  level.  In  the  extreme  case  when  the 
noise  greatly  dominates  the  signal  x^  tends  to  the  trivial 
solutioh:  Xj.-K3. 


Therefore  this  analysis  supports  the  claim  that  the  eigenvectors  are 
perturbed  less  by  noise  than  the  inhomogeneous  solution.  We  point  out  that 
the  validity  of  this  analysis  critically . depends  on  the  validity  of  the 
assumption  that  the  results  for  expected  or  average  cases  are  representative 
of  results  for  actual  cases.  However,  in  the  following  numerical  example, 
actual  cases  are  presented  that  support  the  conclusion  that  the  eigenvectors 
are  perturbed  less  by  noise. 
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2.7  NUMERICAL  EXAMPLE 


Figure  2.1  depicts  the  waveform  that  is  used  in  this  numerical 
example.  The  poles,  ,  and  residues,  Aj ,  which  define  the  uncorrupted 
waveform  are  given  in  Table  2.1.  The  true  number  of  resonances  is  k  ■  12. 
The  time  step  is  At  ■  0.2.  The  waveform  consists  of  400  samples. 


Table  2.2  displays  the  parameter  vector  obtained  with  the  inhomogeneous 
solution  for  the  case  when  n«k  and  M-387  and  for  various  noise  levels.  When 
the  noise  becomes  large  all  parameters  except  c*n  tend  toward  zero.  Table 
2.3  displays  the  weakest  eigenvector,  xNQ£,  for  N-13  and  M*387  and  for  the 
same  noise  levels.  The  parameters  in  this  case  show  no  tendency  toward 
zero.  This  example,  then,  supports  the  conclusion  that  the  weakest  eigen¬ 
vector  is  effected  less  by  noise  than  the  parameter  vector  constructed 
from  the  inhomogeneous  solution. 

Appendix  I  presents  a  technique  to  discriminate  between  the  true  poles 
and  the  extraneous  poles  for  the  case  when  n>k.  This  technique  makes  use  of 
the  tendency  of  extraneous  poles  to  remain  in  the  left  half  of  the  s-plane 
when  the  waveform  is  reversed  while  the  true  poles  are  negated  or  "mirrored" 
through  the  imaginary  axis  of  the  s-plane.  This  tendency  of  the  extraneous 
poles  is  closely  related  to  the  coefficients  tending  to  zero  in  this  example. 
It  follows  that  this  behavior  is  not  observed  if  the  weakest  eigenvector  is 
used  to  estimate  the  resonances.  Appendix  I  describes  in  detail  the  reasons 
behind  this  behavior. 


e.oa 


Table  2.1. 


with 

are 


ices  for  uncorrupted  waveform 
isociated  residues.  Resonances 
iressed  in  terms  of  s-plane  poles. 


IMAG. 

REAL 

0.926 

0.5 

-0.926 

0.5 

2.874 

0.5 

-2.874 

0.5 

4.835 

0.5 

-4.835 

0.5 

6.800 

0.5 

-6.800 

o.: 

8.767 

0.5 

-8.767 

0.5 

10.733 

0.5 

10.733 


0.5 


Table  2.2  Coefficients  for  Prony's  difference 
equation  that  are  constructed  from 
the  inhomogeneous  solution. 


12.56 

-9.724 

6.294 

-3.185 

1.000 


1.442 

1.251 

0.6396 

0.3044 

-1.355 

2.196 

-2.628 

2.450 

-1.914 

1.000 


0.0384 

-0. 0211 


.1452 
.0868 
0.0743 
0.0257 
-0.1543 
0.2438 
-0.2537 
0.1276 
0.0177 
-0.3176 
0.4166 
-0.7654 
1.000 


-0.0112 
0.0007 
-0.0393 
0.0611 
-0.0245 
-0.0607 
-0. 0007 
0. 0192 
-0.0865 
-0.0861 
1.000 
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2.8  INTERPRETATION  OF  NONITERATIVE  ESTIMATION  METHODS  IN  TERMS  OF 

THE  EVD/SVD  ANALYSIS  OF  THE  DATA  MATRIX 

Now  that  the  sense  in  which  the  EVD/SVD  analysis  is  superior  in 
resonance  extraction  has  been  defined,  we  proceed  to  soften  this  claim 
of  superiority  by  demonstrating  that  all  methods  examined  thus  far  can 
be  Interpreted  as  methods  for  constructing  the  parameter  vector  as  a 
weighted  combination  of  the  eigenvectors  of  Q  Q.  The  particular  mode  of 
combination  and  the  particular  choice  of  weights  uniquely  determine  each 
method.  Hence  all  methods  in  the  end  can  be  interpreted  as  particular 
methods  in  which  EVD/SVD  analysis  plays  an  integral  part.  The  question 
of  superiority  reduces  to  the  question  of  what  is  the  best  way  to  combine 
the  eigenvectors  to  form  a  parameter  vector.  Although  we  may  not  come 
to  a  conclusive  answer  to  this  question  in  this  report,  we  examine  several 
first  approaches  for  the  beat  method  of  combination  in  the  section  on  order 
selection  procedures. 

Here  we  will  show  that  our  previous  claim  of  superiority  of  EVD/SVD 
analysis  only  holds  in  a  rather  restricted  case  and  restricted  sense,  and 
that  an  unqualified  claim  of  superiority  in  the  general  case  is  not  wise. 


The  inhomogeneous  solution  can  be  expressed  as  a  linear  combination 
T 

of  the  eigenvectors  of  Q  Q.  To  show  this,  we  begin  by  noting  that  by 
Cramer's  rule  it  is  possible  to  write 
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th  th  T 

where  v^  =  Che  i  column  of  V  =  the  i  eigenvector  of  Q  Q  corresponding 

to  A^.  A  scaled  parameter  vector  constructed  from  x^  which  yields  identical 

pole  estimates  is 


*1  "  [A1N,A2N,,**’ANN]T 

V> 

i.e.,  the  N  column  of  adj  Q  Q.  From  (7)  it  follows  that  x^  is  a  linear 
combination  of  the  eigenvectors: 


N 

Sc, 

i-1  1 


T 

where  the  weighting  coefficients  are  ■  det  Q  Q  v^/A^^  and  v^N  is  the 
element  of  the  i1"*1  row  and  the  N*"*1  column  of  V.  Therefore,  the  estimates 
obtained  with  the  inhomogeneous  solution  are  equivalent  to  the  estimates 
obtained  using  a  weighted  combination  of  the  eigenvectors.  Note  that  the 
eigenvectors  are  weighted  in  inverse  proportion  to  the  corresponding 
eigenvalues. 


In  the  so-called  EVD/SVD  method  of  obtaining  the  parameter  vector 

the  weighting  coefficients  are  all  zero  except  the  coefficient  corresponding 

T 

to  the  weakest  eigenvector  of  Q  Q  which  is  non-zero.  This  method  chooses 
the  weakest  eigenvector  as  the  "natural"  approximation  to  the  eigenvector 
corresponding  to  the  zero  eigenvalue  for  the  noise-free  case. 


The  sensitivity  of:  x^  to  noise  level  is  a  particular  function  of  the 

way  the  eigenvectors  are  weighted  to  form  x^..  For  the  noise-free  case, 

all  coefficients  except  the  one  for  the  weakest  eigenvalue  are  zero  and 

T 

x^  is  simply  a  scaled  version  of  the  weakest  eigenvector  of  Q  Q.  However, 
for  the  noisy  case  x^.  is  not  determined  only  from  the  weakest  eigenvector 
but  is  a  combination  of  ail  eigenvectors.  Since  the  eigenvectors  span  the 
space  of  all  possible  parameter  vectors,  x^  then  has  components  that  may 
not  be  appropriate  for  resonance  estimation.  These  inappropriate  perturb¬ 
ation  of  x  are  maifested  in  the  expected  sensitivity  of  x^  to  the  noise 
level.  We  stress  that  the  components  that  are  appropriate  for  resonance 
extraction  have  no  clear  definition  at  this  time.  Any  absolute  claim  of 
superiority  of  one  method  over  another  is,  in  our  opinion,  ill-advised. 
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So  far  we  have  only  considered  the  case  when  n  -  k  and  W  W  has  only 

one  zero  eigenvalue.  In  this  case,  it  seems  quite  appropriate  to  choose 

T 

the  weakest  eigenvector  of  Q  Q  as  the  "natural"  approximation  to  the  weak- 

T 

est  eigenvector  of  W  W  which  is  the  true  parameter  vector  in  the  noise-free 
case.  But  what  is  the  appropriate  choice  when  k  <  n?  In  this  case,  W  W 

has  more  than  one  zero  eigenvalue  so  that  it  seems  appropriate  to  combine, 
in  some  way,  the  weakest  n  -  k  eigenvectors  of  Q?Q  to  form  a  parameter 
vector.  But  what  is  the  best  method  of  combination?  Also,  it  is  necessary 
to  reduce  the  dimensionality  of  the  parameter  vector  when  k  <  n  since  there 
are  only  k  resonances  and  an  N-dimenaional  parameter  vector  yields  n  ■  N-l 
resonance  estimates.  How  can  this  be  done?  These  questions  will  be  explored 
further  in  the  section  on  order  selection  procedures. 

2.9  JAIN'S  METHOD 

Another  interesting  method  to  find  a  parameter  vector  was  developed 

by  Jain  and  Gupta  [17].  Although  Jain  uses  this  method  with  the  pencil- 
of-functionB  procedure  [18],  we  examine  this  method  in  the  context  of 
Prony's  difference  equation.  The  pencil-of-functions  method  is  described 
in  Appendix  D.  Jain's  method  consists  of  choosing  the  elements  of  x  as 
the  square  roots  of  the  diagonal  elements  of  the  adjoint  matrix  of  Q^-Q. 

The  choice  of  signs  for  the  square  roots  was  not  a  problem  in  the  original 
pencil-of-functions  method  since  for  a  large  class  of  problems  the  square 
roots  were  known  a  priori  to  be  positive.  However,  we  are  not  as  fortunate 
with  Prony's  difference  equation:  there  is  no  means  to  establish  the  signs 
of  the  square  roots  a  priori.  Nevertheless,  it  is  instructive  to  overlook 
this  difficulty  and  to  examine  how  the  method  may  be  interpreted  in  terms 
of  the  EVD  of  Q^Q.  Jain's  method  can  be  interpreted  as  a  means  of  construct¬ 
ing  a  parameter  vector  from  the  eigenvectors  of  QTQ  just  as  the  previous 
two  methods  were  interpreted. 


Jain's  method  chooses  the  parameter  vector  as 


x 


J 
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where  the  choice  of  sign  is  somehow  known  a  priori.  From  (7)  the  diagonal 
T 

elements  of  adj  Q  Q  can  be  expressed  as 

,  det  Q^Q 

Aii  “  2  ~x  id 

j-i  Aj 


where  v^  denotes  the  element  of  the  i*"*1  row  and  column  of  V.  It  is 
clear  then  that  xT  is  a  combination  of  all  the  eigenvectors  even  though 
it  is  not  a  linear  combination.  Again,  each  eigenvector  is  weighted  in 
inverse  proportion  to  the  corresponding  eigenvalue  as  was  done  in  the 
combination  for  x_.  It  follows  that  x  then  is  influenced  by  all  the 

*  J 

eigenvectors  which  may  not  be  appropriate  for  resonance  extraction.  What 

effect  the  mode  of  combination  has  on  x  is  difficult  to  estimate. 

J 


Because  of  the  weights  used  in  the  combination  for  x^,  we  would  expect 
x.j  to  yield  comparable  resonance  estimates  to  those  yielded  by  x^  although 
these  estimates  would  be  different. 


2 . 10  ITERATIVE  TECHNIQUES 

One  can  use  EVD/SVD  of  Q^Q  to  construct  parameter  vectors  and  estimate 
the  resonances.  But  to  obtain  a  fit  to  the  response,  coefficients  A'  must 
be  chosen  to  minimize  (5).  This  two-step  procedure  is  peculiar.  It  would 
seem  that  both  the  z  '  and  A'  should  be  adjusted  simultaneously  for  the 
best  fit  to  the  waveform.  Henderson  [15]  has  shown  that  for  the  noisy  case 
the  resonance  estimates  obtained  by  choosing  the  parameter  vector  as  the 
weakest  eigenvector  of  Q^Q  do  not  provide  the  absolute  minimum  of  (5)  or, 
in  other  words,  they  do  not  provide  an  optimal  fit  to  the  data  in  sense 
of  (5).  The  basis  of  this  behavior  can  be  traced  to  the  fact  that  the 
weakest  eigenvector  minimizes  NQE.  Minimizing  NQE  is  not  equivalent  to 
minimizing  (5). 


To  the  authors’  knowledge,  only  iterative  procedures,  such  as  that 
described  by  Evans  and  Fischl  [19],  (also  described  in  Appendix  G)  can 

Hence, 


be  used  to  choose  both  A'  and  zj,  simultaneously  to  minimize  (5). 

if  we  wish  to  "improve"  the  estimates  beyond  those  obtained  with  the  EVD/SVD 
T 

of  Q  Q  we  must  resort  to  an  iterative  technique.  The  question  of  how  much 


the  choice  of  error  criterion  alters  the  estimates  can  be  approached  with 
a  relation  between  the  two  criteria  obtained  by  Evans  and  Fischl  [19]: 


T 

d  -  X  r 


(8) 


where 


X  - 


a 

n 


0 


•  • 


0  0 


a  i 
n-1 

a 

n  . 


T 

is  an  (M  +  N  -  1  x  M) -dimensional  matrix  and  r  ■  [r. ,r, , . . . ,r.  ]  is 

0  1  nri-n 

the  residual  vector  obtained  by  adjusting  only  the  A  ^  '  to  minimize  (5) 
for  a  given  set  of  z '  corresponding  to  a  given  parameter  vector: 


n 

r .  ■  q.  ”  2  A  (z  )  ,  i  *  0,1,. . . ,m  +  n  . 

j-0  2  2 
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From  (8)  we  observe  that  minimising  the  norm  of  r  provides  a  minimum  upper 
bound  on  the  norm  of  d  through 


T  T 

where  j  | X  j  |  Is  defined  as  the  largest  singular  value  of  X  .  Since  the 

characteristics  of  the  X  matrix  depend  on  the  data,  we  expect  that  certain 
extreme  cases  can  be  found  where  the  two  minimizations  yield  radically 
different  estimates.  Our  own  experience  with  actual  data  indicates  that  the 
two  procedures  yield  comparable  results  when  care  is  taken  to  apply  each 
technique  in  the  proper  fashion.  It  is  clear,  however,  that  minimizing  one 
error  loosely  minimizes  the  other  so  that  if  we  are  not  too  critical  we  may 
claim  that  both  errors  are  approximately  minimized  simultaneously.  Unfor¬ 
tunately  it  is  not  possible  to  predict  how  loosely  the  two  minimizations  are 
coupled  without  knowledge  of  the  particular  data  given  for  resonance  extraction. 


Most  iterative  techniques  attempt  to  minimize  (5) .  Another  novel 
iterative  method,  known  as  the  iterative  generalized  least-squares  procedure, 
does  not  attempt  to  minimize  (5)  but  attempts  to  remove  the  source  of 
"asympotic  bias"  of  the  parameter  vector  by  "whitening  the  residuals".  On 
pages  214-219  of  Eykhoff  [20]  the  source  of  the  bias  is  defined  as 
correlated  residuals.  On  pages  244-247  the  generalized  least-squares 
method  is  defined  and  an  algorithm  due  to  Clarke  [21]  is  presented  as  a 


means  of  implementing  the  method.  Clarke’s  algorithm  is  specialized  and 
applied  to  "whiten  the  residuals"  of  the  least-squares  Prony  technique  in 
Appendix  A  of  this  report.  The  generalized  least-squares  procedure  is 
fundamentally  attempting  to  minimize  the  equation  error  and  is  not  attempting 
to  minimize  (5).  Furthermore,  the  procedure  attempts  to  neutralize  asymtotic 
bias  which  is  only  defined  for  an  infinite  data  sequence  in  the  strictest 
interpretation.  For  these  reasons,  the  generalized  least-squares  procedure 
is  fundamentally  different  from  the  other  type  of  iterative  methods  wa  have 


For  the  inhomogeneous  solution,  which  Clarke's  procedure  begins  with, 

_T- 

the  bias  is  due  the  presence  of  a  second-order  noise  term,  E  E,  in 

QTQ  -  (W  +  E)T(W  +  E)  -  WTW  +  WTE  +  ETW  +  ETE.  This  same  term  is  the  source 
2 

of  the  Mo  I  term  in  (6)  which  was  held  to  cause  the  noise  level  sensitivity 
of  x  .  Perhaps  the  most  significant  comparison  that  can  be  made  is  between 
the  generalized  least-squares  procedure  and  choosing  the  weakest  eigenvector 
as  the  parameter  vector.  Both  of  these  estimation  techqniques  neutralize 
the  ill  effects  caused  by  the  second-order  noise  term  mentioned  above. 
Although  they  will,  in  general,  yield  different  estimates,  the  estimates  are 
expected  to  be  comparable  in  accuracy. 

One  major  drawback  of  the  iterative  procedures  is  that  they  often  do 
not  converge.  Guaranteed  convergence  is  a  requirement  for  any  automatic 
resonance  extraction  procedure.  We,  therefore,  recommend  that  the  iterative 
techniques  be  avoided  in  this  application.  The  use  of  EVD/SVD  analysis  to 
form  resonance  estimates  is  the  safest  procedure  and  can  achieve  the  degree 
of  accuracy  of  the  iterative  techniques. 

2.11  GENERALIZATIONS 

Up  to  this  point,  we  have  considered  only  the  case  that  satisfies 
the  following  restrictions: 

1.  Prony’s  difference  equation  is  used  to  estimate  the  resonances. 

2.  The  impulse  response  of  the  system  can  be  measured. 

New  the  more  general  case  where  these  restrictions  have  been  relaxed  is 
considered.  The  generalization  is  accomplished  by  replacing  Prony’s 
difference  equation  with  the  equation 

n 

iC  “  aj  q^)  -  g±»  1*0,1, . . .  ,m.  (9) 
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where  p  and  are  che  iC^  samples  of  the  filters  on  the  input  and 
output,  respectively,  of  the  system  as  shown  in  Figure  2.2. 

Figure  2.2  depicts  a  single-input,  single-output  linear  system  whose 
transfer  function  is  denoted  HC  and  is  referred  to  as  the  "true  transfer 
function".  The  input  to  the  system  is  denoted  pC  and  is  not  necessarily 
an  impulse.  The  output  to  the  system  is  denoted  qC.  Both  the  input  and 
output  are  corrupted  due  to  measurement  error  to  form  the  measured  input, 
pm,  and  the  measured  output,  qm.  An  estimate  of  the  true  transfer  function 
is  to  be  constructed.  The  poles  and  residues  at  the  estimated  transfer 
function  are  estimates  of  the  resonances  and  residues  of  (1)  for  the  case 
of  nonimpulse  excitation. 


Equation  (9)  can  be  written  in  matrix  form  as  ft  9  *  g  where 


q00  qio 

q01  qll 


q0m  qim 


qn0  P00 
qnl  P01 


q  P 
nm  Om 


nO 

?nl 


nm 


is  the  (M  x  2N) -dimensional  data  matrix  and  0  =  [-a  ,-a. , . . .  ,-a  , 0rt,8. , . .  .3  ] 

u  i  n  u  1  n 

is  the  2N-dimensional  parameter  vector. 

The  filters  in  the  model  of  Figure  2.2  are  typically  chosen  as  first- 
order  discrete  filters  whose  transfer  functions,  F  (z) ,  may  possess  one 
zero  or  one  pole  or  both.  Once  0  is  adjusted  to  minimize  g  an  estimate  of 
the  transfer  function  of  the  system  can  be  formed  as 


H 


m 


8n  Fn  +  ...  +  S  F 
0  0 _ n  n 

aA  F  +  . . .  +  a  F 
0  0  n  n 
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Estimates  of  the  z-plane  poles  of  the  system,  z '  ,  may  be  taken  as  the 
poles  of  H°.  The  Af  may  be  taken  as  the  residues  of  Hm. 


The  special  case  of  Prony's  method  results  from  this  general  scheme 
when  the  F^(z)»z^  and  the  input  is  zero,  i.e.,  p^-0  for  1-0,..., n.  In 
this  case,  Q6-g  reduces  to  Qx-d. 


Any  of  the  methods  for  constucting  a  parameter  vector  with  Prony’s 
difference  equation  can  be  applied,  with  suitable  modifications,  to  construct 
0.  The  most  straightforward  way  to  construct  0  is  to  choose  it  to  minimize 
g  g/9  0  which  is  equivalent  to  choosing  0  as  the  weakest  eigenvector  of  ft  *  ft 
where  the  asterisk  denotes  the  transpose  conjugate. 


■  'jti-  Z 
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2.12  SUMMARY  AND  CONCLUSIONS 


Some  of  the  more  significant  conclusions  of  this  section  are  summarized 
below: 

1.  The  analysis  of  the  noise  effects  on  the  inhomogeneous  solution 

T 

and  the  lowest  eigenvector  of  Q  Q  indicates  that  the  lowest 
eigenvector  is  less  sensitive  to  noise  than  the  inhomogeneous 
solution.  Hence,  the  inhomogeneous  solution  is  considered  less 
desirable  of  the  two  for  resonance  estimation. 

i 

2.  All  noniterative  methods  for  finding  a  parameter  vector;  that  is, 

the  Inhomogeneous  solution,  Jain's  method,  and  £VD/SVD  analysis; 

can  be  interpreted  as  methods  of  combining  the  eigenvectors  of 
T 

Q  Q  to  construct  a  parameter  vector. 

3.  Iterative  techniques  are  fundamentally  different  from  the  nonitera¬ 
tive  techniques  which  we  have  examined  in  this  section.  Most 
Iterative  techniques  choose  the  parameter  vector  to  minimize  the 
"true  error"  between  the  given  data  and  an  approximation  to  the 
data  while  noniterative  techniques  attempt  to  minimize  "equation 
error".  There  is  no  reason  to  believe  that  minimizing  "equation 
error"  produces  less  accurate  estimates  than  minimizing  "true 
error".  While  it  may  be  more  satisfying  to  minimize  "true  error", 
it  is  questionable  that  the  added  expense  and  possibility  of 
divergence  associated  with  an  iterative  technique  is  merited  by  any 
gain  in  accuracy  that  might  exist. 

4.  The  iterative  generalized  lease-squares  algorithm  due  to  Clarke 
[21]  is  fundamentally  different  from  the  other  iterative  techniques 
reviewed  in  this  section  in  that  it  does  not  attempt  to  minimize 
"true  error".  It  is  not  completely  clear  what  this  method  accomplishes 
since  the  criterion  it  attempts  to  satisfy,  that  is,  "white  residuals", 
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cannot  be  realized  exactly  with  a  finite  data  sequence.  If  we  were 
to  use  an  iterative  technique  for  our  application  we  would  prefer 
that  it  satisfy  some  type  of  tangible  criterion  s>«=h  as  minimum 
"true  error".  For  this  reason  the  iterative  generalized  least- 
squares  procedure  is  judged  to  be  less  useful  that  „he  other 
iterative  techniques  for  our  application. 

All  methods  for  constructing  a  parameter  vector  for  Prony's 
difference  equation  carry  over  to  the  generalized  mode.'  with 
suitable  modifications  where  necessary. 

For  the  purpose  of  automatic  resonance  extraction  when  it  is 
known  that  n-k,  the  estimation  procedure  where  the  parameter  vector 
is  chosen  as  the  weakest  eigenvector  of  QTQ  or  n  *  n  i8  judged  to 
be  best  of  all  iterative  and  noniterative  methods.  A  similar  pro¬ 
cedure  for  the  case  where  n>k  is  defined  in  the  section  on  model 
order  selection  procedures. 


3.0  PROCEDURES  FOR  ORDER  SELECTION 


3 . 1  INTRODUCTION 

In  Section  2,  we  assumed  that  k,  the  true  number  of  resonances, 
which  Is  equivalent  to  the  rank  of  W,  was  known.  In  practical  applications 
of  resonance  extraction  procedures  this  assumption  is  unrealistic.  The 
estimation  of  k  for  noisy  data  is  not  a  simple  problem  since  Q  will  be  full 
rank  or  of  rank  N,  whereas  if  there  were  no  noise,  Q  ■  W,  and  the  rank  of 
W  is  easily  determined. 

Furthermore,  in  Section  2,  we  assumed  that  n  ■  k.  In  this  section 
we  examine  the  case  where  n>  k.  Since  the  parameter  vector  is  of  length 
N  »  n  +  1,  there  will  be  n  resonance  estimates  produced  by  the  estimation 
techniques  of  Section  2.  There  will  then  be  n  -  k  extraneous  resonance 
estimates.  We  will  present  a  procedure  developed  by  Henderson  [15]  for 
solving  this  problem  in  this  section. 

In  Section  2,  it  is  shown  that  the  noniterative  procedures  for 

T 

estimating  the  resonances  can  be  interpreted  in  terms  of  the  EVD  of  Q  Q. 
Each  estimation  method  that  we  studied  combined  the  eigenvectors  in  a 
particular  way  to  form  a  parameter  vector.  It  is  not  surprising  that  EVD 
analysis  holds  a  central  position  with  regard  to  all  of  these  particular 
methods  since  EVD  is  simply  a  way  of  rearranging  the  information  in  a 
matrix  in  a  form  that  reveals  the  essential  characteristics  of  the  infor¬ 
mation.  Each  estimation  method  is  a  particular  way  of  combining  the 
information  to  form  resonance  estimates.  For  this  reason  we  adopt  EVD 
analysis  as  a  standard  tool  to  be  applied  in  any  noniterative  automatic 
resonance  extraction  procedure.  Iterative  procedures  are  not  considered 
for  this  application  because  of  their  inherent  problems  such  as  lack  of 
convergence  and  computational  expense. 
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When  n  ■  k  we  concluded  in  Section  2  that  the  weakest  eigenvector  of 

T  T 

Q  Q  was  the  most  "natural"  approximation  to  the  null  eigenvector  of  W  W 

where  the  term  "null  eigenvector"  denotes  an  eigenvector  corresponding  to 

a  zero  eigenvalue.  But  in  the  case  where  n  >  It,  there  are  n  -  k  null 
T 

eigenvectors  of  W  W.  Are  we  then  to  take  some  combination  of  the  n  -  k 

T 

weakest  eigenvectors  of  Q  Q  as  the  most  "natural"  approximation  to  the 
"true"  parameter  vector?  This  question  is  addressed  in  this  section. 

Distributed  linear  systems  such  as  a  scatterer  often  have  an  infinite 
number  of  modes  which  means  that  k  is  infinite.  What  courses  of  action 
are  available  to  us  in  this  case?  In  practical  cases  there  will  be  only  a 
finite  number  of  strong  modes  present  in  the  data  due  to  the  fact  that  in 
any  practical  case  the  dominant  frequencies  present  in  the  excitation  will 
be  limited  to  a  frequency  band  of  finite  width.  Highly  damped  modes  will 
also  be  weak.  The  weaker  modes  will  be  lost  in  the  measurement  noise  so 
that  even  though  k  may  be,  in  truth,  infinite,  practical  considerations 
dictate  that  k  be  estimated  at  a  finite  value. 

This  section  concentrates  on  answering  three  questions: 

1.  How  can  the  true  number  of  resonances  be  estimated  when  n  >  k, 
k  is  unknown,  and  the  noise  level  is  unknown?  (Order  selection 
procedures) 

T 

2.  How  can  the  eigenvectors  of  Q  Q  be  best  combined  to  form  a 
parameter  vector  once  an  estimate  of  the  true  order  is 
available?  (Procedures  for  constructing  resonance  estimates) 

3.  How  can  the  extraneous  resonance  estimates  be  eliminated? 

Only  the  simple  case  where  Prony’s  difference  equation  can  be  applied 
is  considered  in  this  section. 


3.2  STATISTICS  OF  THE  EIGENVALUES  OF  QTQ 

T 

The  problem  of  estimating  the  order  or  rank  of  W  W  reduces  to  the 
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problem  of  discriminating  which  eigenvalues  of  Q  Q  correspond  to  zero 

T  T 

eigenvalues  of  W  W.  If  we  assume  that  M  _>  N  then  Q  Q  will  have  no  zero 

eigenvalues  due  to  the  noise  in  the  data.  In  order  to  isolate  the 

eigenvalues  that  would  normally  be  zero  we  must  come  to  some  understanding 

of  how  the  noise  perturbs  the  eigenvalues.  The  first  step  toward  this 
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understanding  is  realizing  that  the  eigenvalues  of  Q  Q  can  be  described 
as  random  variables.  The  statistics  of  these  random  variables  depends 
to  a  large  degree  on  the  statistics  of  the  noice  in  the  data.  In  most 
practical  circumstances  the  statistics  of  the  noise  will  remain  unknown. 
Beyond  this,  it  is  simply  convenient  to  assume,  in  order  to  remain  completely 
general,  that  the  noie  is  altogether  unknown  but  that  it  has  certain  very 
general  characteristics  by  which  it  can  be  distinguished  as  noise.  These 
characteristics  are  taken  to  be: 

1.  Uncorrelated  from  sample  to  sample  or  white.  (If  we 
assume  that  the  noise  is  due  to  a  very  large  or  infinite 
number  of  very  small  uncuutrolled  influences  which  are 
sufficiently  distant  from  each  other  to  remain  independent 
then  the  noise  can  be  expected  to  be  white) . 

2.  Normally-distributed  with  zero  mean  and  standard  deviation 
of  o.  (By  the  Central  Limit  Theorem  [22]  if  a  random 
variable  is  the  sum  of  an  infinite  number  of  independent, 
zero-mean  random  variables,  each  with  any  distribution 
whatsoever,  subject  to  certain  very  general  constraints, 
then  that  random  variable  is  normally-distributed  and 
zero-mean.) 

The  assumption  of  white  noise  is  believed  to  be  the  best  assumption 
that  can  be  made  when  nothing  is  known  about  the  noise.  This  conclusion 


is  derived  from  the  following  line  of  logic:  If  the  noise  is  correlated 
or  nonwhite,  it  can  be  decomposed  into  a  deterministic  component  and  a 
purely  random  or  uncorrelated  component.  The  deterministic  component 
of  the  correlated  noise  represents  undesired  information  which  cannot  be 
distinguished  from  the  desired  information.  The  deterministic  component 
must  then  be  modeled  by  increasing  the  order  of  the  difference  equation 
model  beyond  that  required  by  the  desired  information.  The  unmodeled 
portion  of  the  noise  is  then  uncorrelated  and  comprises  the  residuals 
of  the  model.  An  uncorrelated  noise  sequence  represents  pure  randomness 
or  lack  of  information  and  is  the  most  difficult  of  all  sequences  to  model 
Summarizing  this  argument,  if  we  are  given  only  the  data  and  absolutely 
no  other  information,  we  must  assume  that  all  deterministic  components 
in  the  data  contain  information  and  the  residual,  purely  random  component 
is  white  noise.  We  desire  to  develop  a  procedure  to  estimate  the  "true" 
order  based  on  the  assumption  that  the  noise  corrupting  the  "true"  data 
is  white. 

Given  these  assumptions  about  the  statistics  of  the  noise,  a 
description  of  the  statistics  of  the  eigenvalues  is  desired.  The  statis¬ 
tical  literature  falls  considerably  short  of  an  adequate  analytical 
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description  of  the  statistics  of  the  eigenvalues  of  a  matrix  such  as  Q  Q 
although  crude  approximations  to  the  statistics  can  be  found  [23]. 
Derivation  of  an  analytic  expression  for  the  statistics  is  not  a  simple 
task.  Our  approach  is  to  derive  an  adequate  approximation  for  the  eigen¬ 
value  statistics  based  on  empirical  studies.  (An  example  of  the  type  of 
studies  performed  is  shown  in  Figure  3.2.) 

Before  defining  a  model  of  the  eigenvalue  statistics  (or  eigen- 
statistics  for  short)  it  is  useful  to  examine  the  eigenvalues  of  the  expec 

tation  or  average  value  of  Q  Q.  If  the  EVD  of  W  W  -  V  A"  V  ,  then  the 

T  W  W  2  WT  N 

EVD  of  £[Q  Q]  ■  V  [A„  +  IMo  ]V  .  In  other  words  the  EVD  of  the  expecta- 
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tion  of  Q  Q  is  equivalent  to  the  EVD  of  W  W  except  that  each  eigenvalue 
2 

has  Mo  added  to  it.  This  is  illustrated  in  Figure  3.1.  The  observed 
eigenvalues  in  specific  random  picks  of  Q  have  some  distribution  about 
these  expected  values.  We  desire  to  find  an  approximation  for  this 
distribution. 
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The  eigenvalue  matrix  has  the  form 


L  -  ^  J 

WWW  T 

where  1-j.  —  ^2  —  *'*  —  are  e^8enva^ues  °f  W  w.  The  model  that  is 

proposed  for  the  eigenstatistics  can  be  written  as  a  particular  form  of 

the  EVD  of  Q^Q: 

T  W  T 

Q  Q  -  V[aJJ  +  Y]  V1 

where  V  is  an  (N  x  N) -dimensional  orthogonal  matrix  whose  columns  are  the 
eigenvectors  of  Q  Q,  Y  is  an  (N  x  N) -dimensional  diagonal  matrix  whose 
diagonal  elements  are  samples  of  the  random  variable  y.  The  frequency 
function  for  y  is  denoted  f(y)  and  completely  describes  the  eigenstatistics 
of  the  model.  This  model  is  chosen  because  it  provides  an  adequate  descrip¬ 
tion  of  the  statistics  that  have  been  observed  in  empirical  studies. 


The  frequency  function,  f(y),  can  be  approximated  as 

2M/N 


fx(y)  -  k 


4By 


(1  +  By)' 


where  K  is  a  scalar  constant  chosen  so  that 


and 


f(y)dy  -  1, 


B  = 


(M  -  3N )0 2 


3-6 


>..•  ».•  ■  in:-., w*v :  j--"'  tiV '1'lV 


■l‘ lV-  sUh£iaii Vi: Wf  V'/ fcmtii ilAduJt liflrt  1  1  r  r~'ii 


The  form  of  this  frequency  function  is  based  on  the  frequency  function 
found  on  page  261  of  reference  [23].  Substantial  modifications  were 
required  to  fit  that  frequency  to  empirically  observed  statistics. 

This  approximation  or  rule  for  the  eigenstatistics  is  quite  accurate 

when  M/N  >  2,  is  acceptable  in  accuracy  when  3/2  <  M/N  <  2,  and  is  generally 

increasingly  unacceptable  in  accuracy  as  M/N  approaches  4/3  or  falls  below 

4/3.  In  practical  situations  quite  often  M  >>  N  so  that  the  rule  usually 

is  sufficiently  accurate.  Figure  3.2  compares  this  rule  against  the 

observed  statistics  for  a  particular  value  of  M/N.  The  histogram  of 

Figure  3.2  is  formed  from  the  results  of  200  Monte  Carlo  trials  where  each 

trial  consists  of  the  EVD  of  Q  Q  and  Q  is  composed  entirely  of  noise 

samples.  The  noise  was  uncorrelated  and  normally-distributed  with  zero- 

mean  and  standard  deviation  a.  All  thirteen  eigenvalues  from  each  trial 

were  used  to  compile  the  histogram  so  that  the  histogram  represents  a 

compilation  of  2600  eigenvalues.  Since  Q  is  composed  entirely  of  noise, 

T  W 

the  eigenvalues  of  W  W  are  all  zero  in  this  case  or  *  0.  It  follows 

T 

that  the  eigenvalues  of  Q  Q  by  our  model  consist  solely  of  samples  of  the 
random  variable  y.  The  approximation  to  the  frequency  function  is  plotted 
with  the  histogram  for  comparison  in  Figure  3.2.  The  rule  has  been  tested 
at  various  other  values  of  M  and  N  with  good  comparisons  in  each  case 
provided  that  M/N  >  3/2.  For  the  case  shown  we  note  that  the  rule  seemr 
to  describe  all  features  of  the  eigenstatistics. 

The  model  of  eigenstatistics  provides  information  about  how  the 

eigenvalues  of  QTQ  should  behave.  In  particular,  the  model  can  be  used 

T 

to  distinguish  those  eigenvalues  of  Q  Q  that  correspond  to  zero  eigenvalues 
of  W^W  since  these  eigenvalues  consist  entirely  of  samples  of  the  random 
variable  y  and  no  other  component,  whereas,  the  other  eigenvalues  have 
deterministic  components  that  serve  to  distinguish  them.  The  method  in 
which  the  model  is  applied  to  discriminate  between  the  eigenvalues  is 
covered  in  our  description  of  order  selection  procedures  which  follows. 


FREQUENCY 


3.3  PROCEDURES  FOR  ESTIMATING  THE  TRUE  ORDER 


Here  we  consider  the  problem  of  estimating  k,  the  "true"  number  of 
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resonances  given  the  EVD  of  Q  Q  for  some  value  of  n  >  k.  Inaccurate 
resonance  estimates  result  when  n  is  less  than  the  number  of  dominant 
resonances  since  in  this  case  the  unmodeled  resonances  appear  as  a  large 
noise  component.  It  must  be  assumed  that  n  can  be  chosen  sufficiently 
large  for  the  resonance  extraction  procedures  to  be  effective. 

In  practical  cases  the  noise  level  is  unknown.  The  order  selection 
procedure  we  use  then  should  not  depend  on  this  information.  We  have  found 
three  procedures  that  are  suited  for  this  problem: 

1.  The  likelihood  ratio  criterion  [24], 

2.  Akaike’s  Information  Criterion  [24,25]. 

3.  Maximum  likelihood  criterion  developed  by  ETI. 

The  first  two  criteria  were  designed  for  use  with  the  minimum  quadratic 

error  for  a  particular  model  order.  Use  of  these  criteria  required  that 

a  model  be  fitted  to  the  data  for  all  orders  l  of  possible  interest. 

Because  fitting  a  model  at  large  number  of  possible  values  for  i  is  an 

expensive  procedure  and  because  we  desire  to  use  the  EVD  to  determine  the 

model  order,  we  adapted  the  first  two  criteria  to  use  the  eigenvalues  of 
T 

Q  Q  instead  of  the  QE.  What  effect  this  modification  has  cn  the  results 
will  be  discussed  when  each  criteria  is  described. 


The  likelihood  ratio  criterion  after  being  adapted  for  the  use  with 
the  eigenvalues,  can  be  stated  simply  as 


1  - 


1  +  1 


<  <S 


This  criterion  is  based  on  the  principle  that  the  eigenvalues 

x 

corresponding  to  zero  eigenvalues  of  W  W  will  be  nearly  equal,  especially 
when  M/N  >■>  1.  6  is  the  threshold  below  which  the  eigenvalues  are  considered 
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equal.  The  order  is  chosen  somewhere  in  the  vicinity  of  where  the 
criterion  is  satisfied  for  several  consecutive  values  of  £ .  The  criterion 
will  in  general  not  continuously  increase  as  Z  increases  because  the 
random  nature  of  the  eigenvalues  is  such  that  the  test  must  be  successful 
several  times  to  be  sure  that  the  proper  order  has  been  reached.  Use 
of  the  criterion  is  illustrated  in  Figure  3.3.  The  order  is  chosen  as 
k'  *  N  -  Z  where  i'  is  the  value  of  Z  for  which  the  criterion  is  satisfied. 
The  6  parameter  could  be  chosen  by  use  of  the  eigenvalues  model  developed 
in  this  section. 

Akaike 's  criterion  as  it  has  been  adapted  for  use  with  the  eigen¬ 
values  is 

2i 

AIC(f)  «  An  XN„£+1  +  m  +  N  -  1 

The  criterion  is  evaluated  for  Z  *  1,...,N.  If  the  value  of  Z  that 
minimized  AIC  is  denoted  Z',  the  estimated  order  is  then  k'  »  N  -Z'  . 

In  order  to  adapt  this  criterion  for  our  purposes,  we  have  probably  taken 
a  few  more  liberties  than  Akaike  would  have  permitted.  If  the  criterion 
defined  above  does  not  work  well  it  should  not  be  taken  as  a  true 
reflection  of  the  performance  of  the  original  criterion.  Figure  3.4 
illustrates  the  application  of  Akaike' s  criterion. 

The  third  order  selection  procedure  is  the  maximum  likelihood  (ML) 

criterion.  This  criterion  was  developed  specifically  for  application  to 

the  eigenvalues  of  Q  Q.  This  procedure  examines  the  smallest  Z  eigenvalues 

and  assesses  the  likelihood  that  those  eigenvalues  are  due  to  noise  only  or, 

T 

in  other  words,  correspond  to  the  zero  eigenvalues  of  W  W.  This  is  done 
by  comparing  the  observed  eigenvalues  to  the  known  frequency  function  of  y. 

The  standard  measure  of  how  well  a  set  of  observed  values  match  a 
known  distribution  is  the  likelihood  function.  The  likelihood  function 
that  is  used  for  this  criterion  is  the  mean  log-likelihood  defined  as 
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Because  the  standard  deviation  of  the  noise  is  unknown,  o  „  must  be 

est 

adjusted  to  maximize  L  .  The  maximum  of  L  is  taken  to  be  the  likelihood 

A*  A* 

that  the  l  smallest  eigenvalues  are  due  to  noise.  The  value  of  o  that 

cS  L 

maximizes  L„  is  an  extimate  of  the  unknown  noise  level.  When  the  maximum 
l 

of  falls  below  a  certain  value,  which  is  called  the  cut-off  likelihood 

and  is  denoted  L  ,  the  observed  eigenvalues  are  so  widely  distributed  that 
c 

the  extreme  values  fall  well  into  the  tails  of  the  frequency  function 

regardless  of  how  a  is  adjusted.  Presumably  when  this  occurs  some  of 

est 

the  eigenvalues  are  not  due  solely  to  noise  and  deviate  significantly  from 
the  noise  eigenvalues. 


The  criterion  can  be  interpreted  in  the  following  way.  a  is 

6SC 

adjusted  to  fit  the  set  of  observed  eigenvalues  into  the  statistical 
distribution  for  noise  eigenvalues.  When  the  best  fit  possible  does  not 
yield  a  sufficiently  high  likelihood  that  the  eigenvalues  are  due  only  to 
noise  then  the  criterion  is  triggered. 
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The  maximum  likelihood  procedure  for  order  selection  consists  of 
the  following  steps: 

1.  For  l  -  2,3,... ,N  adjust  a  ^  to  maximize  L. . 

2.  Choose  l '  as  the  largest  value  of  i  such  that 

max 

°est  Li  *  Lc  for  1  ~  2*3*-*-»£* 

3.  The  selected  order  is  then  N  -  l'  and  the  value  of  a 

est 

that  maximizes  is  an  estimate  of  the  noise  level, o  . 

4.  If  L  >  L  for  l  ■  2,...,N,  then  the  procedure  cannot 

GS  t  jL  C 

discern  a  significant  difference  between  the  observed  eigen¬ 
values  and  the  eigenvalues  of  a  matrix  with  purely  random 
data. 

The  choice  of  the  cut-off  likelihood,  L  ,  determines  how  sure  one 

c 

wants  to  be  that  the  lower  eigenvalues  are  due  to  something  other  than 
noise.  When  the  observed  eigenvalues  are  due  to  noise  only,  the  maximum 
of  tends  to  vary  randomly  about  zero.  Figure  3.5  shows  a  histogram  of 
the  maximum  of  L«.  From  the  histogram  it  can  be  concluded  that  if  L  is 

C 

chosen  as  -1  then  there  is  approximately  one  chance  in  a  hundred  that  the 

criterion  will  be  triggered  due  to  noise  only.  Because  of  the  nature  of 

the  likelihood  function,  this  false  alarm  rate  for  L  -  -1  should  be 

c 

relatively  indifferent  to  the  values  of  M,  N,  or  l.  Lc  *  -1  should 
therefore  represent  a  false  alarm  rate  of  approximately  one  in  a  hundred 
for  all  problems. 

3.4  PROCEDURES  FOR  CONSTRUCTING  A  PARAMETER  VECTOR 

Now  that  several  possible  procedures  for  selecting  the  order  have 
been  defined,  we  consider  the  questions: 

T 

1.  How  can  the  eigenvectors  of  Q  Q  be  best  combined  to  form  a 
parameter  vector  for  resonance  estimation? 

2.  How  can  the  extraneous  resonance  estimates  be  eliminated? 
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In  Section  2,  it  is  demonstrated  that  all  the  noniterative  resonance 
estimation  techniques,  e.g.,  the  inhomogeneous  solution,  Jain's  method, 
or  choosing  the  weakest  eigenvector,  can  be  interpreted  as  different  ways 
of  combining  the  eigenvectors  of  Q  Q  to  form  resonance  estimates.  We 
concluded  that  choosing  the  weakest  eigenvector  as  the  "best"  estimate  of 
the  null  eigenvector  of  W  W  was  the  most  "natural"  estimation  procedure 
for  the  case  when  n  ■  k.  For  the  case  when  n  >  k  the  "natural"  estimate 
of  the  parameter  vector  is  less  clearly  defined. 


Any  vector,  x  ,  in  the  space  spanned  by  the  l  weakest  eigenvectors 


of  Q  Q  satisfies 

I 

tr,  ,rs 

x  >  ii&JlifixJ. 

l  -  ,  T,  a, 

(x  )  (x  ) 

where  N  -  l  is  the  estimate  of  the  "ture"  order.  Stated  differently,  any 

vector  in  the  space  spanned  by  the  l  weakest  eigenvectors  produces  a  value 

£ 

of  NQE  less  than  the  value  of  \  .  It  follows  that  any  x  can  be  consiuered 

AJ 

as  an  approximate  solution  to  Prony'a  difference  equation.  Unfortunately, 

£  £ 

there  are  many  choices  for  x  ,  In  addition,  x  is  N-dimensional  so  that 
it  will  yield  n  resonances  when  N  -  Z  is  the  appropriate  number  of  reson¬ 
ances.  Henderson  [15]  devised  a  method  to  resolve  these  difficulties. 

fch 

Henderson's  method  is  based  on  a  concept  he  termed  "the  l  auto¬ 
regression  nullspace"  which  is  defined  as  the  rowspace  (the  space  spanned 
by  the  rows)  of  the  (&  x  N) -dimensional  matrix 


'«0  al 
n  a0 


0 

a. 


0 

0 


0  a0  at! 


V 


where  the  are,  in  this  case,  the  elements  of  the  true  parameter  vector, 
the  true  order  is  kv  and  N  -  k  +  Z.  Let  h  -  [h^h^.h^, . . .  .h^]1  be  any 
vector  in  the  row  space  of  G  .  Henderson  proved  that  the  roots  of  the 
polynomial , 
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(10) 


constitute  a  superset  of  the  roots  of  the  "true"  polynomial, 


E  Y  "  °- 

3-0 

The  roots  of  the  true  polynomial  are  the  exact  2-plane  poles  of  the  system. 
The  other  2,-1  roots  of  (10)  represent  extraneous  roots.  A  vector  h  can  be 
found  to  produce  any  specified  set  of  extraneous  roots. 

The  essence  of  Henderson's  method  consists  of  approximating  the 

T 

rowspace  of  by  the  space  spanned  by  the  l'  weakest  eigenvectors  of  Q  Q. 
By  making  this  approximation  it  is  possible  to  form  an  estimate  of  the 
"reduced"  (k"*+l) -dimensional  parameter  vector.  Here  N  ■  k'+l*  where  k” 
is  an  estimate  of  the  true  order  k.  The  procedure  is: 

1.  Form  a  matrix  whose  rows  consist  of  the  V  weakest  eigenvectors 
of  QTQ. 

2.  Use  a  Gaussian  elimination  procedure  on  the  matrix  of  step  1 
to  form  an  (2,"*  x  N) -dimensional  matrix  with  zeros  in  the 
opposing  corners.  The  matrix  has  the  form: 


XXXX0000 

OXXXXOOO 

ooxxxxoo 

oooxxxxo 

ooooxxxx 


where  "X"  denotes  a  nonzero  element.  This  matrix  is  intended 
to  approximate  G^  in  some  sense. 

Form  an  (2/  x  k'  +  1) -dimensional  matrix,  H  ,  by  eliminating  the 
zeros  and  appropriately  shifting  the  rows  of  the  matrix  in  step  2. 
Each  row  of  this  matrix  should  be  an  approximation  to  the  true 
parameter  vector  in  some  sense. 
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This  procedure  is  an  ingenious  way  of  getting  around  the  difficulties 
we  mentioned  previously  and  is  intuitively  pleasing.  However,  if  we  choose 
to  be  critical,  it  is  clear  that  some  very  perplexing  questions  could  be 
asked  about  this  procedure.  For  example,  what  error  criterion  or  other 
quality  measure  of  the  estimates  does  the  reduced  parameter  vector  obtained 
with  this  procedure  satisfy? 

On  further  study  of  this  procedure,  other  unanswered  questions  arise: 

1.  The  procedure  begins  by  loading  the  eigenvectors  into  a  matrix. 
Perhaps  the  eigenvectors  should  be  weighted  somehow  according 
to  their  relative  merits  in  estimating  the  resonances  before 
being  loaded  into  the  matrix.  The  weights  might  be  based  on 
how  well  the  corresponding  eigenvalues  fit  within  distribution 
for  the  noise  eigenvalues  in  the  third  step  of  the  ML  proce¬ 
dure  or  on  the  values  of  the  corresponding  eigenvalues  them¬ 
selves.  How  should  the  weights  be  defined? 


v 


2.  The  use  of  Gaussian  elimination  to  form  an  estimate  of 

may  not  be  the  best  procedure  to  use  since  Gaussian  elimination 
can  emphasize  certain  directions  in  the  space  spanned  by  the 
eigenvectors  that  were  not  actually  dominant  originally.  This 
emphasis  of  certain  directions  could  unduly  perturb  the  reduced 
parameter  vector.  It  may  be  possible  to  use  certain  orthogonal 
transformations  in  place  of  Gaussian  elimination  so  that  the 
magnitudes  associated  with  the  directions  are  not  altered.  What 
specific  procedures  can  be  applied  in  place  of  Gaussian  elimination? 

These  questions  merit  further  study. 


3 . 5  NUMERICAL  EXAMPLES 

The  waveform  used  in  the  numerical  example  of  Section  2  is  used 

again  for  these  numerical  examples  which  illustrate  the  order  selection 

procedures  and  the  procedure  to  eliminate  extraneous  resonances.  The 

waveform  is  corrupted  with  uncorrelated,  normally-distributed  noise  with 

zero-mean  and  standard  deviation  o  ■  0.01.  The  data  were  used  to  fill  the 

(M  x  N) -dimensional  matrix  Q  where  M  and  N  were  chosen  to  be  379  and  20, 

respectively.  The  number  of  samples  used  in  filling  QisM  +  N-1*  398. 

The  last  two  samples  of  the  waveform  were  unused.  The  eigenvalues  resulting 
T 

from  the  EVD  of  Q  Q  are  shown  in  Figure  3.6.  We  note  that  It  is  readily 
apparent  to  the  eye  which  eigenvalues  are  due  solely  to  noise. 

Table  3.1  displays  the  resonance  estimates  obtained  using  the  weakest 
T 

eigenvector  of  Q  Q  and  by  using  Henderson's  procedure  to  eliminate  the 
extraneous  poles.  Since  the  poles  must  occur  in  conjugate  pairs  when  they 
are  complex,  we  only  show  the  poles  in  the  upper  half  of  the  s-plane  (posi¬ 
tive-imaginary  region  of  s-plane) .  The  extraneous  resonances  that  result 
from  the  weakest  eigenvector  are  not  shown.  Henderson's  algorithm  does, 
in  fact,  eliminate  the  extraneous  resonances  and  does  yield  relatively 
accurate  estimates  of  the  true  poles.  But  when  the  estimates  using  Henderson's 
procedure  are  compared  to  those  using  the  weakest  eigenvector,  we  note  that 
Henderson's  procedure  produces  less  accurate  estimates.  The  reasons  behind 
the  less  accurate  results  are  thought  to  be  related  to  some  of  the  unanswered 
questions  about  Henderson's  procedure  which  are  enumerated  in  this  section. 

Table  3.2  displays  the  results  of  applying  the  maximum  likelihood 
order  selection  procedure  to  the  observed  eigenvalues  in  Figure  3.6. 

We  recall  that  this  procedure  selects  an  integer,  l',  as  the  largest 
value  of  £  such  that 

max 

a  LJ  <  L  for  i  =  2,3 , . . .1  . 
est  i  c 
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The  column  in  the  table  labeled  max  L  tabulates  the  maximum  of  L  obtained 

X»  jC 

by  adjusting  o  for  each  value  of  Z.  If  the  cut-off  likelihood  were 
est 

chosen  as  Lc  ■  -1  (which  is  the  recommended  value) ,  then  th»*  procedure 
selects  Z'  ■  8  since  this  value  satisfies  the  above  requirements.  The 
selected  order  is  then  N  -  Z'  ■  12  which  is  exactly  the  true  order.  We 
conclude  that  this  procedure  works  very  well  for  this  particular  example. 

A  more  severe  test  can  be  imagined,  however,  where  the  two  groups  of 
eigenvalues  are  not  so  clearly  distinguishable.  Also,  note  that  the  method 
yields  a  very  reasonable  estimate  of  the  noise  level  at  £  =  8. 

t 

An  attempt  was  made  to  perform  a  more  severe  test  but  representative 
results  were  not  produced.  The  crude  algorithm  we  are  currently  using  to 
search  for  the  maximum  likelihood  breaks  down  for  the  more  severe  case. 

The  problem,  though,  is  purely  mechanical  and  does  not  indicate  a  true  limit¬ 
ation  of  the  method.  The  true  test  o*‘  the  procedure  is  its  application  to 
real,  measured  data. 

Table  3.3  shows  the  results  of  applying  Akaike's  criterion  to  the 
observed  eigenvalues.  The  minimum  value  of  AIC  (Z)  defines  the  order. 

From  Table  3.3,  we  observe  that  the  minimum  occurs  at  l  m  1.  One  wonders 
if  a  larger  value  of  N  were  use.  whether  AIC(£)  vrould  go  through  a  clearly 
defined  minimum.  If  we  assume  that  l'  *  1  is  an  accurate  prediction,  then 
the  selected  order  is  N  -  l '  *  19  which  is  not  the  true  order.  We  conclude 
that  Akaike's  criterion  does  not  work  very  well  when  applied  to  the  eigen¬ 
values.  Thi3  conclusion  does  not  indicate  that  Akaike's  original  criterion 
is  not  effective  since  the  criterion  we  have  used  has  been  considerably 
modified  to  enable  application  to  the  eigenvalues.  Akaike's  original 
criterion  was  never  intended  to  be  used  in  this  way. 

Table  3.4  displays  the  results  of  applying  the  likelihood  ratio 
criterion  to  t^e  observed  eigenvalues.  We  observe  that  the  values  of 
1  ••  for  Z  <  8  do  not  differ  significantly  from  the  values  for 

Z  >  8.  There  is  then  no  clean  break  point  for  choosing  the  order.  We 
conclude  that  the  likelihood  ratio  criterion,  as  ue  have  applied  it,  is 
not  suitable  for  selecting  the  order  from  the  observed  eigenvalues. 
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3.6  SUMMARY  AND  CONCLUSIONS 


Some  of  the  more  significant  results  of  this  section  are  summarized 

below. 

1.  Henderson's  procedure  for  constructing  a  parameter  of  reduced 
dimension  has  been  tested  on  synthetically-generated  data  and 
has  produced  accurate  resonance  estimates  although  the  estimates 
were  less  accurate  than  those  obtained  using  the  weakest  eigen¬ 
vector  procedure.  Some  unanswered  questions  about  Henderson's 
procedure  merit  further  study. 

2.  Three  methods  for  estimating  the  true  number  of  resonances  were 
introduced  in  this  section.  On  testing  each  procedure  with 
synthetically-generated  data  we  found  that  only  one  of  the 
procedures  gave  a  reasonable  estimate  of  the  true  number  of 
resonances.  This  method  was  the  maximum  likelihood  procedure. 

3.  A  model  has  been  developed  for  the  statistics  of  the  eigenvalues 

T 

of  Q  Q  assuming  that  the  noise  corrupting  the  waveform  is  white 
and  gaussian-distributed.  The  model  is  used  in  the  ML  procedure 
for  selecting  the  order. 

4.  The  HFTI  algorithm  [16]  is  judged  to  be  less  useful  for  order 
determination  because  it  provides  no  means  of  constructing 

a  parameter  vector. 

5.  Another  method  for  order  estimation  is  the  time-reversal 
technique  described  in  Appendix  I.  Almost  nothing  is  known 
about  the  precise  way  in  which  this  procedure  is  able  to 
discern  which  resonances  are  extraneous.  We  are,  therefore, 
hesitant  to  recommend  this  procedure  until  it  is  better 
understood. 
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4.0  AUTOMATIC  PROCESSING  CONSIDERATIONS 

’  4.1  THE  AUTOMATIC  RESONANCE  EXTRACTION  PROBLEM  FOR  RADAR  TARGET 

IDENTIFICATION 

There  are  two  types  of  resonance  extraction  problems  that  can  be 
(  associated  with  the  radar  target  identification  application: 

1.  Resonance  extraction  for  target  recognition. 

2.  Resonance  extraction  for  target  intelligence. 

i  All  that  is  required  in  the  first  problem  is  to  select  a  pole  set, 

from  a  library  of  resonances  for  known  targets,  that  best  characterizes 
an  observed  waveform.  The  observed  waveform,  in  this  case,  is  the  time- 
domain  radar  return  from  the  target  of  interest.  This  problem,  depending 
I  on  how  it  is  approached,  can  be  much  simpler  than  the  second  problem.  The 

second  problem  involves  the  identification  and  characterization  of  a 
previously  unknown  target.  In  contrast  to  the  first  problem,  this  problem 
has  no  apriori  information  to  work  from. 

The  first  requirement  in  each  of  these  problems  is  that  the  radar 
return  be  measured  in  some  fashion.  We  discuss  some  preferred  ways  in 
which  the  measurements  should  be  performed  next. 

4.2  THE  DATA  ACQUISITION  PROBLEM 

There  are  three  parameters  that  describe  the  mode  of  measurement 
^  of  a  sampled  transient  waveform  which  are  important  for  the  effectiveness 

of  resonance  extraction  procedures.  These  are  the  period  of  observation, 

T,  the  sampling  rate,  m,  and  the  number  of  samples  measured,  N^.  The 
period  of  observation  is  simply  the  time  duration  of  the  waveform  from 
the  first  sample  to  the  last.  The  sampling  rate  is  defined  by  u  *  2n/At 

> 

where  At  is  the  time  interval  between  successive  samples.  These  parameters 
are  related  by  T  «  2irNu/u). 

M 
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There  are  four  more  parameters  that  describe  a  "region  of  feasibility" 
and  a  "preferred  locus"  in  to  -  T  space  shown  in  Figure  4.1.  These  parameters 


are: 


The  expected  period  of  time  that  the  transient  returns  are 


above  the  noise,  T 


max , 


One-half  cycle  of  the  lowest  frequency  of  interest,  T 


tt/oo. 


min  "'"LOW, 

The  highest  sampling  rate  available  or  feasible  on  the  measurement 

apparatus,  u  ,  and 
’  max 

The  Nyquist  sampling  rate  for  the  highest  frequency  of  interest, 

2u.. 


min 


HIGH* 


If  adequate  measurements  of  the  frequencies  of  interest  are  to  be 
made,  id  and  T  must  fall  within  the  region  of  feasibility.  For  adequate 
measurements  to  be  possible  at  all,  it  is  required  that  the  region  of 
feasibility  exist  or  that  conuay  >  mnWn  .  Also,  the  possibility  of 
adequate  measurements  requires  that  the  maximum  number  of  samples  that 
is  feasible,  (N  )  ,  satisfy 


U!  .  T  , 
min  min 

2tt 


The  preferred  locus  defines  the  most  desirable  portions  in  the  region 
of  feasibility  (and  outside  this  region  as  well).  It  should  be  noticed 
that  the  locus  is  a  broken  line.  The  corners  of  the  locus  have 
special  meanings.  The  corner  at  and  Tmav  represents  the  point 

where  holding  the  sampling  rate  constant  while  increasing  the  number  of 
samples  ceases  to  provide  any  more  information.  Any  further  increase  in 
the  period  of  observation  improves  nothing  since  the  added  portions  of 
the  waveform  are  dominated  by  noise.  Further  information  can  be  provided 
only  by  increasing  the  sampling  rate  while  increasing  the  number  of  samples. 
Thus  the  corner  is  formed.  By  continuing  to  increase  the  number  of  samples 
the  second  corner  is  reached.  This  corner  represents  the  point  where 
further  increase  in  the  sampling  rate  is  impossible.  At  this  point  the 
only  option  available  is  to  increase  the  period  of  observation  (even 
though  this  provides  little  benefit)  while  increasing  the  number  of  samples. 
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The  primary  objective  in  data  acquisition  is  compiling  the  maximum 

amount  of  information  about  the  object  of  interest.  This  can  be  done  best 

if  the  largest  number  of  samples  that  is  feasible  to  measure  and  store  on 

the  measurement  apparatus  is  used.  Therefore,  the  point  on  the  preferred 

locus  for  which  Nw  ■  (N  )  defines  the  best  values  of  ai  and  T  for 
M  M  max 

resonance  extraction. 

We  realize  that  there  is  much  more  :o  the  data  acquisition  problem 
than  simply  selecting  the.  sampling  rate  and  period  of  observation  but 
beyond  just  making  accurate  measurements  these  are  the  only  parameters 
of  concern  for  resonance  extraction 


4.3  TWO  FACTORS  EFFECTING  THE  ACCURACY  OF  PRONY'S  METHOD 


Quite  often  in  the  literature  the  sensitivity  of  Prony's  method  to 
noise  in  the  data  has  been  noted.  The  sources  of  this  sensitivity  can  be 
isolated  to  two  factors:  dense  sampling  and  bias.  The  term  "dense  sampling" 
means  that  the  sampling  rate  of  the  waveform  is  much  higher  than  the 
sampling  rate  required  to  faithfully  record  the  highest  frequency  of 
interest  or  the  Nyquist  frequency.  The  dense  sampling  problem  can  be 
ameliorated  by  proper  manipulation  of  the  data.  The  bias  problem  can  also 
be  ameliorated  and  is  discussed  after  the  dense  sampling  problem. 

4.3.1  THE  DENSE  SAMPLING  PROBLEM 


Kulp  [14]  has  performed  an  excellent  study  of  the  effects  of  the 
sampling  rate  on  the  accuracy  of  the  estimates  obtained  by  Prony's  method. 
The  effect  which  Kulp  demonstrated  can  be  explained  by  the  following  simple 


example.  Consider  the  noise-free  exponential  waveform,  w  *  exp(-at).  Two 
samples  of  this  waveform,  which  are  spaced  At  apart  in  time,  are  given  by 
wq  -  exp  (-at)  and  w1  -  exp (-at  -aAt) .  The  "resonance"  for  this  waveform 
is  an  s-plane  pole  at  s  ■  -a.  The  two  noise  corrupted  samples  q^  and  q^ 
which  correspond  to  w^  and  w^  can  be  used  to  form  an  estimate,  s^,  of  the 
single  resonance  as 


This  expression  is  nothing  more  than  a  first-order  Prony  procedure  to 
estimate  the  single  resonance.  It  can  be  verified  that 


The  point  to  be  made  is  that  the  error  in  the  estimate,  s^  +  a,  is 
inversely  proportional  to  the  sampling  interval 


s^  +  a 


_1_ 

At 


£n 


Vo 
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This  same  point  can  be  made  for  higher  orders  and  for  the  least-squares 
version  of  Prony* s  method  as  well  as  all  procedures  that  are  based  on  Prony' s 
difference  equation.  Another  quantity  that  displays  a  dependence  on  the 
sampling  interval  is  the  condition  number  of  the  matrix  W  W  (in  the  notation 
of  this  report)  which  is  used  in  the  least-squares  version  of  the  Prony 
procedure.  In  fact,  Kulp  used  the  condition  number  to  derive  bounds  on  the 
error  in  the  coefficients  of  Prony 's  difference  equation.  The  dependence 
of  the  condition  number  and  the  dependence  of  the  error  in  this  example 
on  sampling  rate  are  different  aspects  of  the  same  phenomenon.  By  this 
example  we  have  briefly  explained  the  problem  of  dense  sampling. 
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When  we  are  given  a  waveform  with  no  other  information  how  do  we 
know  whether  it  is  densely  sampled  or  not?  One  possible  procedure  is  to 
take  the  Fourier  transform  of  the  waveform  to  determine  its  frequency 
content.  If  most  of  the  energy  is  concentrated  in  the  lower  frequencies, 
the  waveform  is  densely  sampled. 

One  way  of  removing  the  effect  of  dense  sampling  is  to  use  a  pre¬ 
processing  technique  such  as  the  redundant-averaging  scheme  described  in 
Appendix  E  of  this  report.  Another,  perhaps  more  ^irect,  way  is  altering 
the  way  in  which  the  data  matrix  is  filled.  For  example,  if  the  Fourier 
transform  indicates  that  the  sampling  interval  can  be  increased  D  times, 
where  D  is  an  integer,  without  aliasing  any  of  the  dominant  frequencies 
of  the  signal,  then  an  alternate  form  of  the  data  matrix  can  be  used  as 
before  to  estimate  the  resonances: 


where  m  >  D  -  1  is  required  if  all  of  the  data  are  to  be  used.  In  order 
to  obtain  the  proper  s-plane  pole  estimates  At  must  be  multiplied  by  D. 

4.3.2  THE  BIAS  PROBLEM 

The  bias  problem  is  very  complex.  There  are  different  definitions  for 
the  term  "bias".  The  definition  that  we  use  is  the  standard  definition  of 
bias  found  in  basic  textbooks  of  statistics  [22],  If  we  denote  the 
estimated  parameter  vector  as  x  and  the  true  parameter  vector  as  x  then  x  is 
said  to  be  an  unbiased  estimate  if 
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where  x  is  an  N-dimensional  random  variable  whose  random  nature  is  due 
to  the  randomness  of  the  noise  in  the  data.  The  unbiased  nature  of  any 
estimate  is  strongly  dependent  on  the  particular  procedure  used  to  form 
the  estimate. 

The  bias  problem  can  be  observed  and  approached  from  many  directions. 
Eykhoff  [20]  describes  the  "asympotic  bias"  in  terms  of  "correlated 
residuals".  The  bias  can  be  described  in  terms  of  subop timal  parameters 
in  the  sense  that  they  do  not  minimize  the  "true  error"  defined  by  (5) 
although  this  approach  to  the  problem  involves  factors  other  than  just  the 
bias  such  as  the  choice  of  error  criterion.  The  fact  that  equation  error 
of  Prony's  difference  equation  is  minimized  by  the  inhomogeneous  solution 
might  lead  one  to  conclude  that  the  particular  choice  of  equation  error  as 
the  quantity  to  be  minimized  causes  the  bias.  On  closer  examination  we 
must  disagree  with  thi3  conclusion.  It  is  not  the  choice  of  error  criterion 
that  causes  the  bias  but  how  the  error  criterion  is  minimized  or,  equi¬ 
valently,  how  the  minimizing  parameter  vector  is  constructed.  Our  obser¬ 
vations  of  Section  2  seem  to  indicate  that  if  the  parameter  vector  is 
chosen  as  the  weakest  eigenvector  of  the  transpose  product  of  the  data 
matrix,  then  the  best  estimates  are  obtained  for  the  case  when  n  ■  k.  The 
symmetrical  perturbations  expected  in  the  eigenvectors  with  uncorrelated 
noise  seem  to  indicate  that  the  expectation  of  the  weakest  eigenvector 
is  equivalent  to  its  "true"  value  without  noise  so  that  the  weakest 
eigenvector  should  be  unbiased.  At  the  present  time,  we  have  no  rigorous 
proof  that  the  weakest  eigenvector  is  unbiased  but  we  feel  that  it  is 
capable  of  proof. 

If  the  noise  is  Gaussian-distributed  and  uncorrelated,  using  an 
iterative  least-squares  technique  to  minimize  "true  error"  produces  an 
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unbiased  estimate.  Unfortunately,  iterative  techniques  are  not  well- 
suited  for  our  application. 

For  the  case  when  n  >  k,  Henderson's  procedure  (Section  3)  should 
effectively  avoid  the  problems  of  bias  based  on  our  current  understanding 
of  the  bias  problem. 

One  method  of  treating  both  the  bias  problem  and  the  dense  sampling 
problem  while  still  using  the  inhomogeneous  solution  is  to  simply  increase 
the  number  of  poles  modeled  or  n.  This  method  is  described  in  Appendix  B. 
Unfortunately,  this  method  suffers  from  the  problem  of  extraneous  poles. 

We  believe  that  the  previous  so-called  problems  associated  with 
Prony's  method  stem  from  the  failure  to  cure  both  the  bias  problem  and 
the  dense  sampling  problem.  Unfortunately,  curing  one  problem  does  not 
cure  the  other  so  that  it  is  very  easy  to  fall  into  one  trap,  if  not  both. 
If  both  of  these  problems  are  treated,  we  have  found  that  Prony's  method 
produces  excellent  estimates  from  data  with  high  noise  levels. 

4.4  THE  TARGET  RECOGNITION  PROBLEM 

The  target  recognition  problem  can  be  broken  into  roughly  four 
subproblems: 

1.  data  acquisition, 

2.  preprocessing  of  the  data, 

3.  data  transformation  or  characterization,  and 

4.  classification. 

We  have  already  discussed  the  data  acquisition  problem.  The  pre¬ 
processing  problem  might  consist  of  choosing  which  portions  of  the  wave¬ 
form  to  use,  filtering  the  waveform,  decimating  the  waveform,  or  other 
similar  procedures.  The  preprocessing  is  intended  to  condition  the  data 


so  that  the  data  transformation  procedure  can  be  applied  efficiently.  The 
data  transformation  problem  can  be  handled  in  different  ways.  One  means 
of  transforming  the  data  may  be  to  apply  a  resonance  extraction  procedure 
to  transform  the  information  into  a  set  of  poles  and  residues.  After  the 
data  are  transformed,  we  must  decide  to  which  class  the  data  belongs.  In 
the  target  recognition  problem  the  classes  are  defined  by  the  pole  sets 
of  known  targets.  The  classification  problem  then  consists  of  choosing 
which  pole  set,  if  any,  in  the  library  of  pole  sets,  best  fits  the  observed 
radar  return. 

The  data  transformation  problem  need  not  be  treated  as  a  resonance 
extraction  problem.  Several  other  approaches  to  this  problem  were  studied 
by  Miller  [26] and  Include: 

1.  Using  a  linear  predictor  from  a  library  of  linear  predictors 
(each  predictor  corresponds  to  a  pole  set  of  a  known  target) 
to  predict  the  next  values  of  the  observed  waveform  and  then 
using  the  mean-squared  error  between  the  observed  waveform 
and  the  predicted  waveform  as  a  measure  of  the  match  [2]. 

2.  Fitting  each  pole  set  in  the  library  to  the  observed  waveform 
by  choosing  residues  to  minimize  the  mean-squared  error 
between  the  modeled  waveform  and  the  observed  waveform- and  then 
using  the  mean-squared  error  as  a  measure  of  the  match. 

3.  Correlating  the  observed  waveform  with  a  library  of  waveforms 
and  using  the  correlation  coefficient  as  a  measure  of  the  match. 

At  the  present  time,  it  appears  that  fitting  each  pole  set  in  the 
library  to  the  observed  waveform  (the  second  approach  enumerated  above) 
is  the  most  reliable  and  easily  implemented  of  all.  The  processing  needed 
for  this  particular  approach  is  small  compared  to  what  might  be  needed 
for  the  application  of  a  resonance  extraction  procedure.  In  addition  an 
optimal  least-squares  fit  to  the  data  can  always  be  obtained  by  using 
linear,  non-iterative  methods.  An  automatic  procedure  using  this  approach 
might  consist  of  the  following  steps: 
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1.  Data  acquisition. 

2.  Discarding  portions  of  the  waveform  that  do  not  contain 
significant  amounts  jf  energy. 

3.  Fitting  each  pole  set  in  the  library  to  the  data. 

4.  Assigning  probability  values,  P  ,  according  to  how  well 
the  i  pole  set  fits  the  observed  data. 

ttl 

is  the  probability  that  the  actual  target  is  the  i  known  target 
in  the  library.  The  rule  by  which  the  probabilities  are  assigned  would 
most  likely  be  developed  by  a  training  process  where  the  target  identifi¬ 
cation  system  observes  variations  in  the  measure  of  fit  of  each  pole  set 
to  a  known  target  at  various  ranges  and  orientations.  The  training  process 
could  continue  even  after  the  system  becomes  operational  by  using  alternate 
means  of  identifying  the  observed  targets.  If  all  the  probabilities  are 
sufficiently  low  then  the  observed  target  is  declared  an  "unknown  target". 

Further  study  is  needed,  however,  before  a  firm  conclusion  can  be 
made  about  which  approach  is  the  best.  In  particular,  attention  should 
be  given  to  which  procedures  are  most  efficient  with  regard  to  the  mode 
in  which  they  are  implemented,  e.g,,  hardware,  software,  parallel  processing, 
or  sequential  processing,  etc.  Such  questions  will  be  studied  further  in 
Phase  II  of  the  current  contract. 

4.5  THE  TARGET  INTELLIGENCE  PROBLEM 

Like  the  recognition  problem,  the  target  intelligence  problem  can 
be  broken  into  subproblems: 

1.  Data  acquisition 

2.  Preprocessing 

3.  Order  selection 

4.  Construction  of  resonance  estimates 

5.  Estimating  physical  features  of  the  target. 


But  unlike  the  recognition  problem,  there  is  no  library  of  information 
to  work  from.  The  data  acquisition  prob1  has  already  been  discussed. 

In  the  following  paragraphs  we  describe  the  most  probable  form  that  an 
automatic  resonance  extraction  system  will  take  based  on  our  present 
understanding  of  th>-  resonance  extraction  problem. 

The  ultimate  aim  of  the  intelligence  problem  is  to  gain  some  kind  of 
useful  information  about  what  the  unknown  target  is.  For  the  purposes  of 
this  report,  we  must  assume  that  thi*  ir.fr-  .ation  is  going  to  be  derived 
from  the  poles.  Then  the  intelligence  problem  considered  in  this  report 
requires  that  poles  be  extracted  from  the  unknown  target’s  radar  return. 

Tills  could  possibJ.e  be  done  in  real  time.  However,  most  likely,  a  wave¬ 
form  would  be  stored  and  processed  at  a  later  time.  Perhaps  the  stored 
T aveft  rm  would  consist  of  an  average  of  all  the  returns  observed  by  the 
target  identification  system  for  the  unknown  target.  Regardless  of  the 
manner  in  which  the  waveform  is  constructed,  .one  automatic  procedure  to 
perform  pole  extraction  is  requires. 

The  preprocessing  step  for  a  Prony-type  resonance  extraction  procedure 
might  consist  cf  the  fol1 owing: 

1.  Discard  portions  of  the  waveform  that  do  not  contain  significant 
amounts  of  energy. 

2.  Fourier-transform  'FFT)  the  waveform  to  determine  if  the  wave¬ 
form  is  densely  sampled. 

3.  Perform  necessary  adjustments  to  correct  the  dense  sampling 
problem  if  it  exists. 

The  otder-seleccion  and  resonance  extiination  subproblems  can  be 
handled  using  the  EVD  of  the  da-a  matrix,  the  ML  procedure  for  selecting 
the  order,  and  Henderson’s  procedure  for  eliminating  the  extraneous 
resonances.  Certain  unanswered  questions  about  Henderson's  procedure 
and  the  ML  procedure  should  be  explored  before  the  exact  form  of  these 
procedures  is  set. 
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The  fifth  subproblem  will  be  studied  in  Phase  II  of  the  present 
contract.  This  subproblem  consists  of  establishing  a  relationship  between 
the  pole  patterns  and  the  physical  characteristics  of  the  target.  At.  the 
present  time  very  little  is  known  about  this  subproblem. 

4.6  SUMMARY  AND  CONCLUSIONS 

Some  of  the  more  significant  results  of  this  selection  are  summarized 

below: 

1.  Guidelines  have  been  established  for  selecting  the  optimal 
sampling  rate  and  period  of  observation  for  measurement  of  a 
transient  waveform  that  is  to  be  used  in  resonance  extraction 
procedures. 

2.  Two  independent  factors  effecting  the  accuracy  of  estimates  of 
Prony-type  procedures  have  been  identified.  These  factors  are 
"dense  sampling"  and  "bias".  The  effects  of  dense  sampling 
can  be  alleviated  by  preprocessing  techniques  or  by  alternate 
means  of  filling  t/  '  data  matrix.  The  bias  problem  can  be 
ameliorated  by  proper  construction  of  the  parameter  vector 
(Henderson's  procedure). 

3.  Two  distinct  problems  fall  within  the  target  identification 
definition:  target  recognition  and  target  intelligence.  Each 
of  these  problems  must  be  treated  in  a  unique  fashion. 

4.  Tentative  forms  for  procedures  and  systems  to  handle  both  the 
target  recognition  problem  and  the  target  intelligence  problem 
have  been  defined.  Phase  II  of  the  current  contract  will  further 
refine  these  tentative  plans. 
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5.0  CONCLUSIONS  AND  FUTURE  EFFORTS 

5.1  SUMMAFY  OF  VOLUME  I 
t 


At  the  end  of  each  section,  except  the  Introduction,  a  summary  was 
provided  for  the  results  of  that  section.  Below  we  summarize  the  con¬ 
clusions  that  will  serve  to  determine  the  course  of  further  efforts. 


1.  Iterative  techniques  are  not  considered  for  the  automatic 
resonance  extraction  application  because  of  their  inherent 
problems  such  as  lack  of  convergence  and  computational  expense. 

£ 

2.  The  various  noniterative  estimation  procedures  can  be  inter¬ 
preted  as  different  methods  of  combining  the  eigenvectors  of 
the  EVD  of  the  data.  There  are  certain  "preferred"  methods 
of  combining  the  eigenvectors  of  which  Henderson's  procedure 
is  an  example  for  the  case  n  >  k.  The  "preferred"  methods 
are  thought  to  produce  unbiased  estimates  of  the  coefficients 
of  Prony's  difference  equation  although  no  proof  is  available. 

3.  The  ML  procedure  for  selecting  the  model  order  is  capable  of 
selecting  the  proper  order  without  any  knowledge  of  the  noise 
level  whatsoever  by  simply  observing  the  eigenvalues  of  the  data. 

> 

4.  Guidelines  have  been  established  for  selecting  the  optimal 
sampling  rate  and  period  of  observation  for  measurement  of  a 
transient  waveform  that  is  to  be  used  in  resonance  extraction 

5  procedures. 


5.  Methods  for  automatically  sensing  the  problem  of  dense  sampling 
and  for  allaying  its  effects  have  been  proposed. 

6.  Tentative  forms  for  procedures  and  systems  to  handle  both  the 
target  recognition  problem  and  the  target  intelligence  problem 
have  been  defined. 
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5.2  CURRENT  STATUS  OF  AND  UNANSWERED  QUESTIONS  ABOUT  RESONANCE 
EXTRACTION  METHODS 


At  the  present  time  it  is  possible  to  construct  automated,  efficient 
procedures  to: 

1.  Estimate  the  number  of  resonances  present  in  the  data,  and 

2.  Construct  estimates  of  the  "true"  resonances  that  in  some 
sense  fit  the  original  data, 

with  no  knowledge  of  the  noifie  given  whatsoever.  These  procedures  should 
be  quite  tolerant  to  noise.  The  procedures  are  based  on  the  eigenvalue 
analysis  of  the  data. 

The  maximum  likelihood  procedure  for  order  selection,  described  in 
Section  3,  will  perform  the  first  step.  Henderson's  procedure,  also  described 
in  Section  3,  will  perform  the  second  step.  Although  these  procedures  perform 
adequately  at  this  time,  there  are  some  unanswered  questions,  that  if 
addressed,  could  create  new  ways  to  further  improve  performance. 


One  unanswered  question  relates  to  the  statistics  of  eigenvalues 
for  alternate  forms  of  the  data  matrix  that  might  be  used  to  allay  the 
effects  of  dense  sampling.  If  the  eigenvalue  statistie.3  do  change 
significantly,  then  another  model  for  the  statistics  must  be  developed  for 
the  alternate  forms  of  the  data  matrix  and  used  in  the  ML  procedure  for 
order  selection.  Another  question  relates  to  the  proof  that  EVD  allows 
unbiased  estimation  of  the  coefficients  of  Prony's  difference  equation. 
Such  a  proof  is  needed  to  place  the  use  of  EVD  techniques  on  firm  ground. 

Other  unanswered  questions  about  Henderson's  procedure  which  were 
posed  in  Section  3  are: 

1.  How  should  the  weakest  eigenvectors  be  weighted  prior  to  the 
use  of  this  procedure? 

2.  What  specific  procedures  can  be  applied  in  place  of  Gaussian 
elimination  in  this  procedure? 


i 

Further  study  along  the  lines  of  these  questions  will  almost  certainly 
yield  further  improvements  in  the  procedures.  Perhaps  the  most  important 
I  unanswered  question  is  how  well  the  methods  for  preprocessing,  order 

selection,  and  resonance  estimation  will  perform  with  real  measured  data. 

5.3  FUTURE  EFFORTS 

Under  Phase  II  of  the  current  contract  the  concepts  and  procedures 
developed  in  Phase  I  will  be  further  developed  and  plans  for  specific  non- 
cooperative  target  recognition  (NCTR)  systems  will  be  further  refined. 

To  accomplish  the  above  the  following  questions  must  be  addressed: 

1.  What  is  the  relationship  between  the  natural  resonances  and 
the  physical  shape  and  dimensions  of  a  target? 

F 

2.  What  forms  can  a  NCTR  system  assume?  Which  form  is  the  best 
under  what  conditions? 

3.  What  modifications  are  required  if  the  procedures  are  to  be 

^  implemented  in  hardware? 

Answers  to  these  questions  will  be  pursued  by  simulating  possible  NCTR 
systems  with  measured  data  in  the  second  phase  of  this  contract. 

► 
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7.0  LIST  OF  SYMBOLS 


Roman  Symbols: 


A 


j 


true  residues 


estimated  residues 


B 


C 


i 


D 


d 


d 


i 


dNQE 

E 


e 


i 


F 


j 


f 


G 


g 


a  constant  that  determines  the  value 
of  y  that  maximizes 

coefficients  of  a  linear  combination 
of  the  eigenvectors 

integer  that  specifies  how  many  times 
the  sampling  interval  can  be  increased 
without  aliasing  any  of  the  frequencies 
of  interest 

M-dimensional  equation  error  vector 
f~h 

i  element  of  the  equation  error  vector 

particular  equation  error  vector  that 
satisfies  Qx^  -  dNQE 

error  matrix  defined  as  Q-W  "  E 

t  h 

measurement  error  of  the  i  sample 

th 

transfer  function  of  the  j  filter  in 
the  generalized  model 

frequency  function  for  y 

approximation  to  the  frequency  function  f 

(2,  x  N)  -dimensional  matrix  whose  rows  span 
the  4C  1  autoregression  nullspace. 

equation  error  vector  for  the  generalized  model 
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i  element  of  the  equation  error  vector 
for  the  generalized  model 

(£“  x  k^  +  1) -dimensional  matrix  used  in 
Henderson's  procedure. 

estimate  of  the  transfer  function  of  the 
single-input,  single-output  system  which 
models  the  scatterer 

any  N-dimensional  vector  in  the  row  space 


j  element  of  h 

an  identity  matrix  of  appropriate  dimension 
general  index 
general  index 

constant  that  determines  the  amplitude  of 
true  number  of  poles 
estimated  number  of  poles 
estimate  of  the  true  order  k 

cut-off  normalized  mean  log-likelihood  used  in 
the  maximum  likelihood  procedure  for  selecting 
order 


normalized  mean  log-likelihood  function 

expected  normalized  mean  log-likelihood  for 
pure  noise 

index  associated  with  the  lower  eigenvalues 
number  of  rows  of  Q,  M  *  m  +  1 
equals  M-l  (defined  for  convenience) 
number  of  columns  of  Q,  N  *  n  +  1 
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number  of  samples  in  the  measured  waveform 
number  of  samples  in  the  analyzed  waveform 
number  of  poles  modeled 

T 

transpose  product  of  Q  or  Q  Q 

the  probability  that  the  actual  radar  target 
is  the  ic  known  target  in  the  target  library 

Any  (N  x  N) -dimensional  singular  matrix 

(N  x  N) -dimensional  singular  approxiraant  to  P 

the  i1"^  sample  of  the  output  of  the  filter 
on  the  excitation  in  the  generalized  model 

(M  x  N) -dimensional  measured  data  matrix 

(^x  n) -dimensional  version  of  Q  or  Q  with  the 
N  column  removed. 

an  alternate  form  of  Q  that  is  useful  in  treating 
the  deii3e-sampling  problem 

the  n-rank  approximant  to  Q 

t*H 

the  N  column  of  Q 
th 

the  iLI  sample  of  measured  response  of  the 
linear  system 

the  i^  sample  of  the  output  of  the  j  ^  filter 
on  the  response  in  the  generalized  model 

the  (M  +  N-l) -dimensional  residual  vector 

the  iC^  element  of  r 

(M  x  N)-dimensional  diagonal  matrix  whose 
diagonal  elements  are  non-negative  and  are 
called  the  singular  values  of  Q 

true  s-plane  poles 
estimated  s-plane  poles 
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s 

s 


s 

T 


T 

max 

T 

imin 

t 

At 

U 


V 


,w 


w 

w 


X 


x 


x 


the  (M  x  N) -dimensional  matrix  constructed 
from  S  by  forcing  the  smallest  singular  value 
to  zero. 


Laplace  -  transform  variable 

period  of  observation  (length  of  waveform 
in  time) 

expected  period  of  time  the  measured  transient 
waveforms  are  above  the  noise 


one-half  cycle  of  the  lowest  frequency  of 

of  interest,  equals  ir/w. 

Low 


time  variable 


time  step 

(M  x  M) -dimensional  orthogonal  matrix  whoj^e 
columns  consist  of  the  eigenvectors  of  QQ 

(N  x  N) -dimensional  orthogonal  matrix  whose 
columns  are  the  eigenvectors  of  Q  Q 

the  (N  x  N) -dimensional  orthogonal  matri^  whose 
columns  consist  of  the  eigenvectors  of  W  w 

fch  T 

the  i  eigenvector  of  Q  Q  corresponding  to 
or  the  i  column  of  v. 

tijj|?  element  of  V  belonging  to  the  i^  row  and 
j  column 

(M  x  N) -dimensional  true  data  matrix 
uncorrupted  version  of  q 
uncorrupted  value  of  q^ 

matrix  constructed  from  the  coefficients  of 
Prony's  difference  equation 

any  N-dimensional  parameter  vector  or  vector 
of  coefficients  at  Prony's  difference  equation 

aijj^  n-dimensional  version  of  x  or  x  with  the 
N  element  removed 


solution  of  the  inhomogeneous  equation, 

-T — T- 

Q  QXj  -  Q  q 
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*NQE 

*NQE 


J 

a 


Y 


y 


z 


scaled  parameter  vector  constructed  from  x 
which  yields  identical  pole  estimates 

expected  inhomogeneous  solution 

T 

weakest  eigenvector  of  Q  Q 

n-dimensional  vector  construc|gd  from  * 
by  scaling  so  that  its  N^element  :Es 

one  and  thenweliminating  the  n  n  element. 

parameter  vector  constructed  with  Jain's 
method 

any  vector  in  the  space  scanned  by  the'  £ 
weakest  eigenvectors  of  Q  Q 

(N  x  N) -dimensional  matrix  whose  diagonal 
elements  consist  of  samples  of  the  random 
variable  y 

a  random  variable  used  to  model  the  random 
component  in  the  eigenvalues  due  to  noise 
in  the  data 

z-transform  variable 


true  z-plane  poles 


estimated  z-plane  poles 
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Greek  Symbols : 


a 

0 


A 


ij 


6 


e 


A 


N 


coefficients  of  Prony's  difference  equation 

th 

scaling  factor  for  the  output  of  the  j  filter 
on  the  excitation  in  the  generalized  model 

th 

element  of  the  i  row  and  j  column  of 
adj  QtQ 

cut-off  parameter  used  in  the  likelihood  ratio 
criterion  for  selecting  the  order 

i 

2N-dimensional  parameter  vector  for  the 
generalized  model 

(N  x  N) -dimensional  diagonal  matrix  whose  ^ 
diagonal  elements  are  the  eigenvalues  of  Q  Q 


A 


M 


ANS 


A 


W 

N 


X 


1 


X. 


w 

t 


c 


est 


a 


(M  x  M) -dimensional  diagonal  matrix  whose  ^ 
diagonal  elements  are  the  eigenvalues  of  QQ 

diagonal  (N  x  N) -dimensional  matrix  constructed 
from  A^j  by  forcing  the  smallest  diagonal  element 
to  zero 

(N  x  N) -dimensional  diagonal  matrix  whose 
diagonal  matrix  whos^  diagonal  elements  are 
the  eigenvalues  of  vJ  w 

th  X 

i.  eigenvalue  of  Q  Q  the  eigenvalues  are 

ordered  X,  <  X0  <  .  . .  <  X.. 

1  -  i.  -  —  N 

th  T 

i  1  eigenvalue  of  W  w,  the  eigenvalues  are 

ordered  A?  <  X^f  <  .  . .  <  X^ 

1  —l  -  —  N 

expectation  operator 

adjustable  parameter  used  to  form  the  maximum 
likelihood  |stimate  of  cr  from  the  lower  eigen¬ 
values  of  Q  Q 

the  standard  deviation  of  the  assumed  Gaussian- 
distributed,  zero-mean,  and  uncorrelated  noise 
corrupting  the  samples  of  the  waveform 
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^awnwffn 


a 


“high 

“low 

“max 


> 


the  (M  x  2N) -dimensional  data  matrix  for  the 
generalized  model 

sampling  rate  (angular  frequency) 

highest  frequency  of  interest 

lowest  frequency  of  interest 

Maximum  feasible  sampling  rate 

Nyquist  frequency,  equals  200^^ 


* 


* 


> 


I 


> 
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Miscellaneous  symbols: 


a| 


V 

x 


adjA 


detA 

InA 


exp  (A) 


denotes  the  euclidean  norm  of  vector  A 
or  matrix  A 

gradient  operator  with  respect  to  the 
vector  x 

denotes  the  adjoint  matrix  of  matrix  A 
or  the  transpose  of  the  matrix  of  cofactors 
of  A 

denotes  the  determinant  of  matrix  A 

denotes  the  natural  logarithm  of  some  real 
or  complex  number  A 

denotes  e,  the  base  of  natural  logarithms, 
raised  to  the  power  A 

denotes  the  inverse  of  square  matrix  A 

denotes  the  transpose  conjugate  of  complex 
matrix  A 


